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ABSTRACT 


The momentum and continuity equations for a minor gas are combined with 
the momentum equation for the major constituents to obtain the time dependent 
continuity equation for the minor species reflecting a wind field in the back- 
ground gas. This equation is used to study the distributions of helium and argon 
at times of low, medium, and high solar activity for a variety of latitudinal- 
seasonal wind cells. For helium, the exospheric return flow at the higher 
thermospheric temperatures dominates the distribution to the extent that much 
larger latitudinal gradients can be maintained during periods of low solar ac- 
tivity than during periods of high activity. By comparison to the exospheric flow, 
the smoothing effect of horizontal diffusion is almost negligible. The latitudinal 
variation of helium observed by satellite mass spectrometers can be reproduced 
by the effect of a wind system of air rising in the summer hemisphere, flowing 
across the equator with speeds on the order of 100 to 200 m/sec, and descending 
in the winter hemisphere. Argon, being heavier than the mean mass in the lower 
thermosphere, reacts oppositely to helium in that it is enhanced in the summer 
hemisphere and depleted in the winter. By using winds which are effective in the 
lower thermosphere, the anomalous vertical helium profiles observed from 
rockets can be reproduced. The time response of the helium density distribution 
following the initiation of a wind field implies the likelihood of a factor of two to 
four density enhancement at night over the daytime values. 
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NOMENCLATURE 


Aj> nm = coefficient defined in Appendix D 
Bg = coefficient defined in II, B. 2 

Bj? nm = coefficient defined in Appendix D 
b = radius to base of exosphere 

C = coefficient defined in Appendix D 

nm 

D = molecular diffusion coefficient 
g = local acceleration of gravity 

k T 

H = = scale height of minor gas 

mg 

k T 

H 1 = = scale height of major gas 

J = coefficient defined in Section II D 

k = Boltzmann's constant 

K = eddy diffusion coefficient 

m = molecular mass of minor gas 

M = molecular mass of major gas 

n = number density of minor gas 

N = number density of major gas 

p = pressure 

P m (6) = Legendre polynomial 

r = radial coordinate 

S = shape factor in exponential temperature profile 
T = temperature 

t = time 
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= flow velocity of minor gas 
= flow velocity of major gas 
= mean molecular speed = (2.55 k T/m) 1/2 

= thermal diffusion factor 
= factor determining wind velocity gradient 
= /3 1 x 10 2 (defined in III. B.l.a) 

= gamma function = ( l - 1) ! for n = integer > 0 

= coefficient defined in II. D. 

= polar angle (colatitude) 

= cosine 9 
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THERMOSPHERIC WIND EFFECTS ON THE GLOBAL DISTRIBUTION 
OF HELIUM IN THE EARTH'S UPPER ATMOSPHERE 

I. INTRODUCTION 

The enhancement of upper atmospheric helium in the winter hemisphere 
has been noted from satellite mass spectrometers (Reber and Nicolet, 1965; 
Reber, et al., 1968) and has been suggested to explain anomalies in satellite 
drag data (Keating and Prior, 1968; Jacchia, 1968). The best mapping of this 
phenomena has come recently from the quadrupole mass spectrometer flown on 
the OGO-6 satellite (Reber, et al., 1971). Figure 1 shows the distribution of 
helium from this measurement taken over half an orbit near sunrise on 7 June 
1969. As the data are taken over a range of altitudes, this parameter is normal- 
ized out by dividing each measured density by the predicted density from the 
appropriate Jacchia model atmosphere (Jacchia, 1965, hereafter referred to as 
J65), 

C _ Measured helium density 
6 N Model helium density 

To the extent that J65 correctly represents the real temperature profile, and the 
atmosphere is in diffusive equilibrium, [He] is the ratio of the actual helium 
density at 120 km to the constant boundary density of the model. The data in 
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GEOGRAPHIC LATITUDE 


Figure 1. Ratio of the measured helium number density to Jacchia 65 model atmosphere 
density as a function of geographic latitude (Reber, et al., 1971). 

Figure 1 shows that this ratio varies by an order of magnitude over the geo- 
graphic latitude range 50*N to 80*S, with a peak near 55*S. It is further shown 
that the location of the peak in the helium density is closely correlated with the 
geomagnetic dipole field, with peak locations falling quite close to the 53° south 
magnetic dipole latitude. 

Vertical profiles of constituent densities obtained from rocket measure- 
ments consistently show departures from diffusive equilibrium profiles for 
helium and occasionally for argon as well. Kasprzak (1969) summarizes seven 
of these flights and attributes the profiles to an upward flux of helium ranging 
from 2.0 x 10 8 cm -2 sec -1 to 2.6 x 10 10 cm -2 sec -1 . He notes that these values 
are consistent with fluxes calculated by McAfee (1967) for lateral transport of 
helium at the base of the exosphere due to diurnal temperature variations. Hart- 
mann, et al. (1968) reports an enhanced helium distribution in the winter 
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thermosphere and decreased argon concentrations. They attribute their results 
to a lowering of the turbopause level below that assumed by the COSPAR Inter- 
national Reference Atmosphere (CIRA, 1965). Reber (1968) interpretes the de- 
viations from diffusive equilibrium profiles in terms of the long diffusion times 
in the lower atmosphere and the combination of this phenomenon with changes in 
time of either the turbopause level or the exospheric temperature. 

In the discussion of helium data from Explorer 17, Reber and Nicolet (1965) 
suggest that the observed latitudinal/seasonal (spring- fall) variation of a fac- 
tor of two could be explained by a seasonally dependent change of 5 km in the 
turbopause altitude. The variation of helium density with turbopause altitude has 
been studied in detail by Kockarts and Nicolet (1962); Kockarts (1971), exploring 
this mechanism further, points out that a factor of 20 variation in the eddy dif- 
fusion coefficient is required to explain the OGO-6 data. However, Colegrove, 
Hanson and Johnson (1965) comment that the molecular oxygen to atomic oxygen 
ratio at 120 km is proportional to the eddy diffusion coefficient. As an oxygen 
variation of this magnitude is not observed, there is an apparent inconsistency 
in the use of this mechanism to explain the entire helium variation. 

Johnson and Gottlieb (1969, 1970) suggest that the source of the winter 
enhancement of helium is a large scale meridional circulation system, with air 
moving from the summer polar regions toward the winter pole. They infer 
a downward flow on the order of 100 cm/sec between 150 and 200 km altitude 
in the winter polar region from the compressional heating required to maintain 
the temperature in this region. Taking into account the upward flux required to 
support exospheric transport (McAfee, 1967) due to the density enhancement, 
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they arrive at a concentration buildup of about a factor of two at the winter pole. 

They state further that there is probably an inverse effect over the summer pole 
so the circulation mechanism would support a pole-to-pole ratio of about four. 

To investigate the effect of winds on minor constituents in more detail, 

Reber, Mayr and Hays (1970) studied the continuity equation for a minor gas 
modified to include the effects of winds. The simplified wind system used in 
their calculation consists of a constant vertical velocity above 200 km and zero 
wind below that altitude, along with a cosine distribution in latitude. They con- 
clude that the global helium distribution can be explained on the basis of upper 
thermospheric winds and that these winds would affect the vertical distribution 
to altitudes below 100 km. 

In the present work these calculations are expanded in several ways to re- 
flect a more realistic physical situation. Basically, the analytical approach in- 
volves combining the momentum and continuity equations for a minor gas (e.g. 
helium) with the continuity equation for the major background gas (in this case, 
the total of 0, 0^ and N 2 ). The result is the three-dimensional time-dependent 
continuity equation for the minor gas, modified from the usual version by the 
addition of motion in the background gas through which the minor gas is diffus- 
ing. The horizontal component of the meridional wind field is expressed in terms 
of the vertical component through the major gas continuity equation, while the analytic 
form of the vertical wind permits a more realistic and easily varied circulation 
cell to be described. The calculation includes the smoothing effect of horizontal 
diffusion at all altitudes, and in addition, the upper boundary condition reflects 
the exospheric transport discussed by McAfee (1967) and Hodges and Johnson 
(1968). 
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The resulting differential equation is integrated numerically, using an 
IBM 360-75 computer. The results, presented in the following chapters, permit 
the effects of vertical wind profile, exospheric temperature, horizontal diffusion 
and exospheric transport to be examined in detail with respect to their influence 
on the horizontal and vertical distribution of helium. In a later chapter, the types 
of wind fields which produce distributions consistent with satellite and rocket 
observations are described. The majority of the study is carried out for the steady 
state situation (primarily to reduce computer time), but the time response follow- 
ing a sudden initiation of a wind field is investigated for a number of typical 
systems. Finally, the behavior of argon (which exhibits an opposite reaction to 
winds compared to helium) is examined from the point of view of (1) comparison 
with measurements and (2) emphasizing the physical processes important in the 
redistribution of gases in the upper atmosphere. 

II. DYNAMIC MODEL FOR A MINOR GAS 

The problem of studying the three dimensional distribution of a minor gas, 
when a motion field is impressed on the major (background) gas, is approached 
by combining the momentum and continuity equations for the minor gas 
with the continuity equation for the major species. The result is the minor gas 
continuity equation, modified from the usual form by the addition of terms re- 
flecting winds in the major species. 

The calculation is simplified considerably by the assumption (discussed in 
detail in IIB) that the wind fields and minor gas distribution can be averaged over 
a day; thus, any longitudinal variations are neglected. A solution to the continuity 
equation is obtained by expanding the minor gas distribution and the wind field 
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in terms of Legendre polynomials and solving for the coefficients of the gas 
distribution for a given wind field. The remainder of this chapter is devoted to 
a discussion of this calculation; detailed derivations are found in the appendix. 


A. Combined Continuity and Momentum Equations 

In spherical coordinates, with no longitudinal variation and no sources or 
sinks, the continuity equation for the minor species (n) becomes 


Bn 3 , » 2nv r 1 3 

Bt + 3r (0 + r + r sin 6 B 6 


(sin d n Vg) = 0, 


where 

n = number density of minor gas, 


( 1 ) 


t = time, 

v = flow velocity of minor gas, 
r = radial coordinate, 

6 = polar angle. 

Similarly, the radial and latitudinal components of the momentum equation 
become 


n [v r - V r ] = -D 


3n n(l+a)3T n v 
3r + T 3r + H 


3 n n 3 T n 

3 r T 3 r 


( 2 ) 


and 


n [v e -V 6 ] = -5 


3n n (1 + a) 3T 
+ T 30 


(3) 


where 


V = flow velocity of background (major) gas 
D = molecular diffusion coefficient, 
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a = thermal diffusion factor. 


T = temperature, 
kT 

H = — = scale height of minor gas, 
mg 

k = Boltzmann's constant, 

m = molecular mass of minor gas, 

g = local acceleration of gravity, 

K = eddy diffusion coefficient, 
kT 

H* ~ Mg ~ sca ^ e height of major gas, 

M = mean molecular mass of major gas. 

The eddy diffusion term, 


Bn n BT n 
B r T B r jj' 


is added to the expression for the radial momentum by considering the flux to 
be composed of diffusive and eddy components, after the development by Cole- 
grove, et al. (1965). (Horizontal eddy diffusion effects are not included in the 
present calculation.) By combining equations (1), (2) and (3) with the continuity 
equation for the background gas one can obtain a form of the minor gas continuity 
equation which contains the effect of motion in the background gas (see Appendix A): 


3n _ _B_ 

B t B r 


2 

+ — 


(d 

3n n (1 + a) 3T n”l ^ 

Bn n BT n 

1 

l 

3r T d r Hj 

B r T B r 

j 


3 n 

. n (1 + 

a) B T n 

: 4. 

+ K Bn 

n BT 

n 

+ 

B r 

T 

Br H 

B r 

+ T Br 

H'J 


v [” ri d N d n"l „ 1 f"n B N B n"| 
r |_n b7 “ bTJ + 9 7 [_N dd ~ T0J 

a (dn n(l+a) BtN 

| Dsm ® + — T — sej 


r 2 sin 6 


(4) 
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where N = background gas number density. The first two terms on the right 
represent the effect of radial diffusive flow in a spherically symmetric atmosphere; 
solution of the equation containing only these two terms, with the diffusive flux, 
n(v - V. ), set equal to zero yields the usual static diffusive-equilibrium vertical 
distribution. The fifth term reflects the smoothing effect of horizontal diffusion. 

The third and fourth terms represent the perturbations introduced by verti- 
cal and horizontal winds, V r and V e , on the minor gas distribution. The physical 
effect of this type of term can be better seen by using the approximations 


and 


1 3N ^ J. 
N 3r u' 


1 3n ~ 1 

n 3 r H 


With these, the vertical wind term becomes 



( 5 ) 


Thus, an upward wir’d will cause a decrease in density for gases whose mass is 
less than the mean mass, and an increase in density for gases whose mass is 
greater than the mean mass; for gases whose mass is close to the mean mass 
in regions where the wind is important (e.g. N 2 ) there is little effect. The op- 
posite sense holds true, of course, for a vertically downward wind. 

One can also study the reaction to a vertical wind in terms of its effect on 
the composition of a cell of air moving with the wind. A cell moving upward in 
the region where mixing is no longer important will maintain the relative com- 
position of air at a lower altitude. For a light gas such as helium, this results 
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in a decrease in the relative number density at higher altitudes. Again, the 
opposite holds true for a downward wind: a cell of air transported to lower alti- 
tudes reflects the relative composition of the higher altitude, resulting in an en- 
hancement of the lighter species and a depletion of the heavier gases. While the 
vertical wind is distorting the vertical profile from a diffusive distribution, 
molecular diffusion is attempting to re-establish this profile. Thus, to be ef- 
fective, the vertical wind speed must be significant relative to the local diffusion 
velocity, v D , where 

D 

V d “ H ‘ 

In this sense, the process may be considered analogous to the competition be- 
tween eddy and molecular diffusion in establishing the transition from the mixed 
to the diffusive atmosphere (the "turbopause"). 

B. Assumptions and Approximations 
1. Longitudinal averaging 

In the derivation of equation 4 one assumption has already been made and 
noted, that of no longitudinal (or diurnal) variation in the quantities of interest. 
This is done on the basis that the phenomenon being studied is an averaged, 
relatively long term, effect (namely a latitudinal-seasonal phenomenon) as op- 
posed to a diurnal effect. The wind system postulated is also a diurnal average 
and is divergence free in the east-west direction; the only requirement is that 
outflow during the day near the summer pole must exceed the inflow during the 
night, while at the winter pole the inflow must exceed the outflow. A wind sys- 
tem of this general nature is discussed in some detail by Johnson and Gottlieb 
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(1970) as being required to explain the relative warmth of the mesosphere and 
thermosphere over the winter pole in the absence of direct solar heating. 
There is also no indication from available data that there exists a persistent 
longitudinal variation in the helium distribution. 

2. Polynomial expansion of minor gas distribution and wind field 
Equation (4) is solved by expressing the latitudinal wind field and minor 
gas (hereafter specified as helium) distribution as an expansion in Legendre 
polynomials: 


and 


V r (r, 0) = 


n (r, 6 ) = 



\ (r) (9) 


"n(OP n (*)- 


( 5 ) 


( 6 ) 


The horizontal wind components are related to the vertical components through 
the major gas continuity equation. 


Bn b 

B t B r 


2NV 

(NV ) + - 

r 


_1 _B_ 

r sin 9 B 9 


(sin 9NV g ) = 0. 


(7) 


For a steady state the latitudinal components become (see Appendix C) 


where 


V* (r. 0) = - 


H Bj(r )Pi‘(0), 
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B^(r) = r 


By? 


+ r 


N 
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( 9 ) 


and 



r(l) 
rtf + 2> 




r(^ = 



1 e X dx = (*£ - 1) ! 


for n = integer > 0- 
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3. Model atmosphere 

The molecular and eddy diffusion coefficients (D and K), the major gas 
number density (N), and atmospheric temperature (T) are assumed to be inde- 
pendent of latitude. Between 80 kilometers and 120 kilometers the atmospheric 
parameters of number densities and temperature are taken in tabular form from 
the CIRA 1965 model atmosphere (CIRA, 1965). (For the numerical solution, it 
was found desirable to modify the temperature profile slightly to eliminate dis- 
continuities in the slope. (See Appendix B.) Above 120 kilometers the analytic 
expressions for the temperature and major constituent density profiles are 
taken from the 1965 Jacchia model atmosphere as modified by Walker (1965). 

The major (background) gas is taken to consist of molecular nitrogen (N 2 ), 
molecular oxygen (0 2 ) and atomic oxygen (O). 

C. Solution of the Minor Gas Continuity Equation 

Using in Equation (4) the expansions from Section IIB, multiplying each 
term by P m (9) sin0 , and integrating over the polar angle from 0 to tt, we find 
the continuity equation for the m th harmonic in the expansion of the helium dis- 
tribution (for details, see Appendix D): 


2 3n m _ 2 _3_ 

2m + 1 3t 2m + 1 dr 


2 s 2 

+ b — 

2m + 1 nm r 
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where 


and 


= D 


Bn n (1 + a) B T n m 

Tr“ +_ ~T "3r H 


+ K 
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O) P m (/x)d/x # 


sin 

* / -i 


a / BP ( M )\ 

— P (/lz) — - (sin 0 — — ) dfj., 

6 30 26 




= cos 0, 


A numerical solution to Equation (10) is obtained using an integration technique 
described by Lindzen and Kuo (1969). Details on the method of solution are given 
in Appendix D. 

D. Boundary Conditions 

The helium density at the lower boundary (80 km) is taken to be 1.989 x 
10 9 cm" 3 from CIRA,1965, and assumed to be independent of latitude. This im- 
plies that there is a sufficiently large reservoir at this altitude to supply or 
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accept the amount transported horizontally in the thermosphere, with no modi- 
fication of the lower boundary density. 

At the upper boundary (500 km, the base of the exosphere) the slope of the 
helium profile is determined from the vertical flux across this level: 


n(v r - V r ) = -D 


Bn n (1 + a) 3T n 
3 r + T 3 r + H 


- K 


Bn n BT n 
3 r T B r H'_ 


( 11 ) 


where all the quantities are evaluated at 500 km. Above this altitude, molecules 
are assumed to be describing ballistic trajectories and returning without ex- 
periencing collisions. This results in a horizontal flow (exospheric transport), 
related to horizontal temperature and density gradients, which has been studied 
extensively by McAfee (1967) and Hodges and Johnson (1968). The expression 
developed by the latter authors is used here to express the vertical helium 
flow in terms of atmospheric properties at 500 km: 


nv r (l + V 2 (n <v> H 2 ), 

where 



( 12 ) 


< v) = mean molecular speed, and 
b = radius to base of exosphere. 

Going through a development similar to that outlined in IIC, and setting 
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Equation (11) reduces to 


-2 J 


L 


n (n + 1) 
2 n + 1 


S n - 


/ n vp Ap 
L I n 4 'Ll 

n J{/ 


+ D 




n n( 1+a >3T n n\ 
+ T 3r + H ) 


= 0 , 


(13) 


for the coefficient of the m th term at the boundary. (Terms reflecting the effect 
of eddy diffusion are dropped, since at 500 km altitude K« D.) 


HI. MINOR GAS RESPONSE TO LARGE SCALE MOTIONS 
IN THE MAJOR SPECIES 
A. Eddy Diffusion Coefficient 

The individual component density in the upper thermosphere and exosphere 
is extremely sensitive to the value of the eddy diffusion coefficient in the region 
of transition from a mixed atmosphere to one controlled by molecular diffusion 
(Lettau, 1951; Kockarts and Nicolet, 1962; Mange, 1961; Colegrove, Johnson 
and Hanson, 1966). In particular, the effect is enhanced for minor species (such 
as helium) whose mass differs greatly from the mean mass of the mixed at- 
mosphere. Also, the sense of the effect for a particular gas depends on the dif- 
ference in mass between the gas and the mixed mean mass: for a heavier gas 
such as argon, an increase in eddy diffusion coefficient will result in an in- 
creased density at higher altitudes, while for helium an increased eddy diffusion 
coefficient will decrease the density. These effects are discussed in detail 
in the references cited above. 

For the study of a minor gas response to winds, it is necessary to include 
a realistic function for the eddy diffusion coefficient. A number of profiles were 
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Figure 2. Diffusion coefficients as a function of altitude: (A) eddy diffusion coefficients 
used in the calculations presented here; (B) constant eddy diffusion coefficient which pro- 
duces samehigh al ti tude hel ium density as (A); Johnson and Gottlieb (1970) eddy diffusion 
coefficient based on thermal considerations; (D) molecular diffusion coefficient for = 

1200. 

tried, from a constant to an approximation of a profile suggested by Johnson 
and Gottlieb (1970) based on thermal considerations, consisting of a 
constant value over an altitude interval with an exponential decrease above and 
below this interval. These two profiles are shown in Figure 2 along with the 
Johnson and Gottlieb profile and the molecular diffusion coefficient (for an exo- 
spheric temperature of 1200®). The profile A falls off from 130 km rather than 
150 km as suggested by Johnson and Gottlieb due to the lower absolute maximum 
value; they state that the decrease should begin about a scale height above the 
altitude where the eddy diffusion and molecular diffusion are comparable. 
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HELIUM DENSITY AT 500 km (cm" 3 ) 


To determine a realistic value for the eddy diffusion coefficient, helium 
densities obtained from the mass spectrometer on OGO-6 are compared against 
several calculated values for each of the two diffusion coefficient profiles 
(Figure 3). The influence of dynamics on the distribution of helium is most 
likely minimal during quiet periods near equinox. Thus, the mass spectrom- 
eter densities are taken from a magnetically quiet period (A p = 5) on 24 Sep- 
tember 1969 at latitudes of +48° and -41° and at 500 km altitude. The exospheric 
temperature, determined from the molecular nitrogen density (measured by the 
mass spectrometer), was 1176°K and 1235°K at +48° and -41° latitudes respec- 
tively. Helium densities at these two locations were 2.73 x 10 6 cm -3 and 



EDDY DIFFUSION COEF. (cm 2 sec _1 ) 

Figure 3. Helium density at 500 km for = 1200° as a function of eddy diffusion coefficient. 
The two curves correspond to the eddy diffusion coefficient profiles shown in Figure 2. 
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2.19 x 10 6 cm -3 . The calculated values are obtained by using the same com- 
puter program used in the dynamic-diffusion calculations, setting the wind equal 
to zero, and using 1200° for the exospheric temperature. It can be seen from 
Figure 3 that either of the two eddy diffusion profiles mentioned can produce sat- 
isfactory agreement with the high altitude data and one need only choose the appro- 
priate constant value; either 1.8 xlO 6 cm 2 sec -1 for the constant profile or 2.5 xlO 6 
cm 2 sec -1 for the Johnson and Gottlieb profile satisfy the data. 

The shape of the vertical profile for helium is not affected by the eddy dif- 
fusion coefficient used in the calculation. This can be seen in Figure 4 where 

helium density profiles for various constant eddy diffusion coefficients are shown 



HELIUM (cm -3 ) 


Figure 4. Helium density as a function of altitude for various constant values of the 
eddy diffusion coefficient and for the eddy diffusion coefficients of Figure 2 (marked 
A and B). The exospheric temperature, , is 1200°. 
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along with the two profiles which best match the high altitude data (marked A and 
B in the figure). Figure 5 shows the expanded display of the same set of calcu- 
lations where now the dependent variable is the ratio of the number density of 
helium to the sum of the major gases (molecular nitrogen, molecular oxygen and 
atomic oxygen). For the remainder of the calculations, the Johnson and Gottlieb 
eddy profile is used, with a maximum value of 2.5 x 10 6 cm 2 sec* 1 , and a scale 
height of 9.1 km for the exponential regions. It is assumed that the eddy diffusion 
coefficient determined during equinox conditions can be applied globally during 
solstace conditions. 



Figure 5. Relative abundance of helium as a function of altitude for various eddy diffusion 
coefficients; the curves A and B refer to the eddy diffusion coefficient profiles of Figure 2. 
T* = 1200°. 
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B. Effect of Latitudinal — Seasonal Circulation 

There are no direct measurements of vertical velocities in the upper meso- 
sphere and thermosphere, nor are there any measurements of large scale lati- 
tudinal-seasonal circulation cells in these regions. That high winds exist in the 
thermosphere and mesosphere is in little doubt, however, and many analytical 
studies have been published concerning various aspects of upper air wind sys- 
tems based on assumed pressure gradients (Geisler, 1966, 1967; Dickinson and 
Geisler, 1968; Dickinson, Lagos, and Newell, 1968; Lindzen, 1966, 1967; Chap- 
man and Lindzen, 1970; Volland, 1969; Volland and Mayr, 1970; Mayr and Volland, 
1971; Kohl and King, 1967a, 1967b, 1968; Challinor, 1968, 1969; Bailey, Moffett 
and Rishbeth, 1969; Rishbeth, Moffett and Bailey, 1969) or thermal require- 
ments (Johnson and Gottlieb, 1970). A number of general features of these studies 
have been extracted and incorporated into a large scale, easily parameterized 
circulation cell: (1) the air flows from regions of relatively high pressure to 
regions of low pressure (in this case, from the summer to the winter hemisphere); 
(2) the vertical profile of the vertical component of the wind field consists of a 
rapid increase with altitude up to heights where viscous effects may be expected 
to become important, at which point the velocity tends toward a constant value 
(the profile and magnitude are consistent with those of Kohl and King (1967) and 
Volland and Mayr (1970)); (3) the horizontal component of the wind field is 
determined from the vertical component through the major gas continuity equa- 
tion (see Sec. IIB and Appendix C). Studying the effect of thermospheric wind 
cells of this type on the distribution of minor species is the object of the present 
work. This investigation comprises three general areas: (1) a broad parametric 
steady state analysis corresponding to periods of low, medium and high solar 
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activity, and demonstrating the effects of such phenomena as exospheric trans- 
port, horizontal diffusion and wind profiles; (2) a specific comparison with 
OGO-6 mass spectrometer helium measurements near the June 1969 solstace, 
showing the types of wind systems which are compatible with the observations; 
(3) the results of time-dependent calculations showing rates of generation of 

various minor gas distributions when a wind system is suddenly "tumed-on". 

1. Cellular Motion 

a. Vertical Profile 

The adoption of an arbitrary, easily parameterized wind system is fa- 
cilitated by the reduction of the horizontal and vertical components of this sys- 
tem to expressions in terms of the vertical components only (see Section IIB and 
Appendix C). This reduction, accomplished by using the major gas continuity 
equation and expanding the wind field in Legendre polynomials, greatly simpli- 
fies the analysis as it permits a complete description of a circulation cell with 
a minimum number of parameters. 

The general vertical velocity profile imposed on the circulation cells con- 
sists of a rapid increase with height up to altitudes where viscous effects begin 
to dominate, at which point the velocity tends to become constant with altitude. 
This shape can be expressed as 

Wp 

Ye =2 (1 +erf^(z -z^)]}, (14) 

where = vertical major gas velocity ^-£5^, 

W ^ = maximum value of > 

erf(x) = error function (erf(x) = 2//7T J q x e -t dt), 
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z 0 l ~ reference altitude (where slope, dv/3z, is equal to 1 //tT w/ 3 or where 
V^W^.km), 

/3£ = factor determining altitude gradient (equal to ^ /wp 3 /dz a t 
z = z^) (km -1 ). 

Thus, by defining W^, z 0 ^ and /3£, the complete circulation cell is determined 
(for a given density profile; see below). The generalized altitude profile is shown 
in Figure 6; Figure 7 shows representative vertical velocity profiles for w = 100 
cm sec" 1 , z 0 - 200 km and several values for ^ (/?£ = /?£ x 10 2 ). Also shown 
for comparison is the vertical profile deduced by Johnson and Gottleib (1970) from 
thermal considerations for the winter mesosphere and thermosphere. 

b. Cell Shape 

Using the assumptions of Section II.B.l, the air motion is approximated by 
a pole-to-pole circulation cell, with air rising in the summer hemisphere 
(0*0* 90°), flowing across the equator and descending in the winter hemisphere 
(90° < d * 180°). This distribution can be described by using only the first 
Legendre polynomial in the expansion of the vertical component of motion (see 
Section H.B.2.): 

V r (r, 0) = Vj (r) cos 9. (15) 


The latitudinal component then becomes (Appendix C). 


V d (r, d) = - 1/2 Bj (r) sin 0, 


where 


B x (r) = r 


3V, 




(16) 


(17) 
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ALTITUDE (km) 


400 



Figure 6. General vertical profile of the vertical wind used in the calculations. This 
specific profile is characterized by T ^ = 1100°, W = 200 cm/sec, Zg =230 km, and 
/6 = 1.8 x 10“ 2 km -1 
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Figure 7. Vertical wind profiles for several values of fi, with =1000°, W = 100 cm/sec 
and Zq =200 km. Also shown for comparison is the vertical wind profile deduced by Johnson 
and Gottlieb (1970) from thermal considerations. 

Thus, the complete circulation cell is defined in terms of cosine and sine func- 
tions and the profiles of the vertical wind component and major gas number 
density. 

It can be seen by inspection of equation (16) that the direction of latitudinal 
flow will be determined by a balance of the density gradient term on the right 
(always negative) against the positive first and third terms. In regions where 
the density gradient term dominates, the flow will be toward increasing values 
of 0(i.e., toward the winter pole); in regions where the wind gradient and ampli- 
tude dominates, V e (r, 6 ) is negative and the flow is toward decreasing 9 (the 
summer pole). Figure 8 shows the vertical profiles of the vertical and horizontal 
components for a typical wind system, with w = 100 cm sec -1 , z Q = 200 km, 

/3 = 2.0 and an exospheric temperature, T ■, of 1100°. (Henceforth, 
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T co = l 100° 
100/200 
£= 2.0 
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VERTICAL 

(cm/sec) 


WIND VELOCITY 

igure 8. Vertical and horizontal wind profi les for =1100°, W = 100 cm/ sec, Zq = 200 km 
nd 13 - 2.0. Note that the horizontal wind becomes negative (toward the summer pole) between 
40 and 185 km. 
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Figure 9. Direction of the wind vectors associated with the vertical 
and horizontal profiles of Figure 8. 





this nomenclature will be abbreviated, i.e. 100/200, / 3 = 2.0). The region of 
return flow, toward decreasing 6 , is seen to lie between 140 km and 185 km; the 
sharp break in horizontal velocity at 120 km is due to the change in slope of the 
density of the atmospheric model used. The shape of the circulation cell associ- 
ated with these profiles is depicted in Figure 9. 

Qualitatively, the relationship between the gradient of the vertical velocity 
component and the direction and amplitude of the horizontal flow can be seen by 
considering the requirements for continuity in a vertical column. If the vertical 
velocity, v, increases with altitude at exactly the same rate at which the density, 
N, decreases, the flux, Nv, is constant along the length of the column and there 
is no horizontal inflow or outflow (assuming that a diffusive vertical profile is 
maintained for the major species). If, however, the vertical velocity increases 
more rapidly than l/N, the flux out the top of a small volume element in the 
column exceeds the flux coming in through the bottom and there is a need for 
compensating inflow through the sides of the volume element. Conversely, if V 
increases less rapidly than l/N there is a net horizontal outflow. 

The relationship of the altitude regime of reverse flow to the vertical pro- 
file parameters is shown in Figures 10, 11 and 12 for exospheric temperatures 
of 800°, 1100° and 1500°. These three temperatures were chosen as they approxi- 
mate global averages for periods of low, medium and high solar activity. In 
these figures, the region of reverse flow is the area to the right of a particular 
P /altitude curve; the area to the left of a given curve indicates flow toward the 
winter pole. In general, as the height increases at which the vertical velocity 
levels off, there is a corresponding increase in the altitude of return flow. The 
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Figure 10. Altitude regions of reverse flow for Zq = 180, 200 and 230 km; =800°. The 
horizontal wind is toward the summer pole for values in the altitude -/S plane to the right of 
a given curve. 
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small variation in shape of these curves with exospheric temperature reflects 
the variation in the major constituent density profile which is balanced against 
the wind gradient in the calculation of Bj (r). 

2. Minor Gas Response to Cellular Motion: Helium 

The response of minor gases to large scale dynamic systems shows up most 
dramatically in two ways that can be compared directly with observations: (1) 
the vertical profile as would be measured from a rocket borne mass spectrom- 
eter; and (2) large scale latitudinal distributions which can be compared with 
satellite observations. This discussion will emphasize helium as there are 
many more data applicable to this gas; argon will be discussed separately in a 
later section. 

As tne majority of the data on the large scale distribution of helium has 
been obtained at satellite altitudes, it is desirable to compare the results of the 
calculation directly to high altitude data. Figure 13 shows the calculated results 
of the helium density at fixed altitudes of 300 km and 500 km as a function of 
latitude for a typical wind system. It is seen that the distributions are 
smooth functions, increasing from the summer pole toward the winter pole, with 
a pole-to-pole ratio (R p ) of 8.7 at 500 km. (This parameter, the pole-to-pole 
density ratio, turns out to be a useful quality figure, and will be referred to 
frequently in later sections.) 

Figure 14 displays the vertical helium number density profiles associated 
with the latitudinal distributions of the previous figure. Three profiles are given, 
representing the summer pole (0°), the winter pole (180*) and for comparison, 
the static-diffusion profile. In general, for the simple wind cells studied, the 
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HELIUM (cm-3) 



Figure 13. Helium density as a function of polar angle for the constant altitudes of 300 km 
and 500 km. The exospheric temperature is 800° (low solar conditions), W = 50 cm/sec, 
Z 0 = 200 km, /3 = 1 8. Also shown are the densities in the absence of winds. 

polar profiles represent the extrema of the dynamic effects on the helium 
density. 

The following three sections will examine the effects of exospheric lateral 
transport, exospheric temperature, and horizontal diffusion on the vertical and 
horizontal distributions of helium. Following these is a discussion of the 
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ALT (km) 



HELIUM (cm' 3 ) 


Figure 14. Helium density at the summer and winter poles as a function of altitude, for the 
50/200, [i- 1.8 wind system of Figure 13. Also shown is the static profile. 

variations in the distribution as functions of the shape, amplitude and altitude of 
the wind cell. 

a. Exospheric Transport 

By far the largest effect tending to smooth out horizontal variations in 
helium density is that due to lateral flow in the exosphere. This transport is 
proportional to the quantity 



(see equations (12) and (13), section n.D) where <v> is the mean molecular speed 
and e = b/H, where b is the geocentric distance and H = kT/mg. It can be seen 
that an increase in exospheric temperature will cause an increase in J t through 
both the mean molecular speed and the scale height, H. J is shown as a function 
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J (cm/sec) 


of exospheric temperature in Figure 15, where it will be observed that an in- 
crease of more than a factor of five in exospheric flow has resulted from a 
temperature increase from 800° to 1500°. The net result of this transport 
mechanism is a flux up and out of the thermosphere in regions of comparatively 
high helium density (i.e. near the winter pole), high altitude flow over the equator 
toward the summer pole, and flow into the thermosphere from the top in the 
summer hemisphere. 



Figure 15. The exospheric transfer velocity function, J, as a function of exospheric 
temperature. The exospheric flux is related to J through the expression shown. 


32 




(lUVi 005 ) d H 


There are no measurements of this exospheric flux, but its existence can 
hardly be in doubt. The magnitude of the transport, however, might be questioned, 
so this quantity was varied from 0 up to 1.5 times the magnitude suggested by 
Hodges and Johnson. The results for two wind systems are shown in Figure 16, 
where it can be seen that the pole-to-pole ratios (R p ) vary only slightly when the 
jxospheric flow is varied in the neighborhood of the Hodges and Johnson value. 



: igure 16. The pole-to-pole helium density ratio, R p , at 500 km as a function of exospheric 
'ux for average solar conditions (Tg, = 1100°). The value 1.00 $ represents the value col- 
lated from Equation 12 for the Hodges and Johnson (1968) flux. 
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When the exospheric flow is removed, however, the pole-to-pole ratios in- 
crease by more than an order of magnitude. The smoothing effect on the large 
scale distribution can be seen in Figure 17, where the latitudinal variations of 
helium at altitudes of 120 km, 300 km and 500 km are shown for; (1) no exo- 
spheric flux and for; (2) half the Hodges and Johnson value. For the "no flux" 
case the higher harmonics are clearly present, while for the other case they are 
gone. The vertical profiles corresponding to these situations are shown in 
Figure 18. The calculations discussed in the remainder of this paper are 
performed using the Hodges and Johnson expression for the magnitude of the 
exospheric transport. 



Figure 18. Helium density vertical profiles at the poles for no exospheric flux and half the 
Hodges and Johnson flux. Shown for comparison is the no wind helium profile. 
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b. Horizontal Diffusion 


Lateral diffusion below the base of the exosphere is the other "restoring 
force" acting to smooth out horizontal variations in component density, but com- 
pared to exospheric flow its effect is minor. Figure 19 shows the variation in 
helium density at 500 km with latitude for a typical wind system, both with and 
without horizontal diffusion included in the calculation. The effect on the pole- 
to-pole ratio is 10%, certainly less than could be observed with present measuring 



Figure 19. Helium density at 500 km versus polar angle for the case of no 
horizontal diffusion and including horizontal diffusion. 
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techniques. At 120 km the pole-to-pole ratio is increased by 16% by eliminating 
horizontal diffusion, but here again it would be difficult to observe. Figure 20 
gives a comparison of vertical profiles with and without back diffusion. The 
amplitude of this effect would grow with increasing atmospheric temperature 
due to the temperature dependence of the molecular diffusion coefficient, but its 
relative importance as a smoothing mechanism would diminish compared to the 
much more temperature sensitive exospheric transport. 


c. Exospheric Temperature 

It is to be expected that increasing the amplitude of the winds in a circula- 
tion cell (as might be anticipated during periods of high solar activity) will produce 
a corresponding increase in their influence on the distribution of minor gases. 
Increasing the atmospheric temperature, however, has the opposite result, as the 



Figure 20. Helium density vertical profiles at the poles corresponding to Figure 19. 

Also shown is the static profile. 
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exospheric transport increases strongly with temperature and tends to ffeduce 
any variation. The competition between these two effects is such that at periods 
of high solar activity, much stronger winds are required to produce and main- 
tain a given latitudinal variation than at periods of low solar activity. This re- 
lationship is illustrated in Figure 21 where the pole-to-pole ratios of helium 
(R p ) at altitudes of 120 km, 300 km and 500 km are shown as functions of the 
maximum vertical velocity in a cell; three exospheric temperatures are repre- 
sented, corresponding again to periods of low, medium and high solar activity. 

A number of interesting features are apparent here. First, a factor of two 
in R p at 500 km requires an order of magnitude higher wind for an exospheric 



w(cm/sec) 


Figure 21. Pole-to-pole ratios, R p , at 120, 300 and 500 km as functions of maximum vertical 
wind speed, w. Low, medium and high solar conditions are represented; Zg = 200 km and 
j3- 1.8 for all the curves. 
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temperature of 1500° than for one of 800°, in agreement with the discussion in 
the previous paragraph. Second, for higher temperatures the latitudinal vari- 
ation is suppressed at high altitude relative to 120 km. This is another conse- 
quence of the enhancement in exospheric return flow at high temperatures, 
which tends to smooth out latitudinal variations most significantly at higher 
altitudes. Thus, at times of low solar activity there should be much better 
agreement between low altitude measurements (e.g. from rockets) and satellite 
measurements than during periods of high solar activity. Also, it is not likely 
that the wind amplitude increases sufficiently (due only to increases in the 
pressure gradients), from periods of low to high solar activity, to maintain the 
same level of disturbance in the helium distribution; therefore, the observed 
pole-to-pole ratio should be highest at times of low solar activity. 

d. Shape, Amplitude and Altitude of Wind Cell 

While the effects discussed so far, particularly the exospheric transport, 
have an important bearing on the ultimate distribution of a minor gas, it is the 
wind field itself which sets up the variation from a uniform, static-diffusion dis- 
tribution. In this section, we shall examine the effect of varying the character- 
istics of the wind field itself, specifically the altitude and amplitude of the cell 
and the altitude of the return flow. 

The pole-to-pole ratios of the helium density as functions of maximum ver- 
tical wind speed are shown in Figures 22 through 27 for /?’ s of 1.8 and 
4.0, z 0 *s of 170 km, 200 km and 230 km, and exospheric temperatures 
of 800°, 1100°, and 1500°. The corresponding vertical profiles at the summer 
and winter poles and equator are given in Figures 28 through 33; Figures 28, 31 
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Figure 22. Pole-to-pole ratios, R , at 120, 300 and 500 km versus vertical wind speed 
w, for T* =800°, Z 0 = 180, 200 and 230 km, and 1.8. 
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Figure 23. Pole-to-pole ratios, R , at 120 300 and 500 km versus vertical wind speed 
w, for T ra =800°, 2 Lq =180, 200 and 230 km, and /3 = 4.0, 

and 33 compare equal vertical amplitudes, while Figures 29, 30 and 32 give the 

profiles for approximately equal values of R p . The latitudinal distributions at 

120 km, 300 km and 500 km are given in Figures 34 through 39 for the same set 

of parameters. Representative horizontal and vertical wind profiles for these 

systems are shown in Figures 40 through 43, and Figures 44 through 47 picture 

the corresponding circulation cells. The plots of R p versus wind speed indicate 
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Figure 24. Pole-to-pole ratios, R , versus vertical wind speed, w, for T_ =1100 and /3 = 1 .8. 
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Figure 25. Pole-to-pole ratios, R , versus vertical wind speed, w, for =1100 and / 3 = 4.0. 
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Figure 30. Helium density versus altitude for T^ -1100° and the wind systems /3 - 1.8, 
70/180, 90/200, 130/230 which produce similar values of R p (500 km). 
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Figure 35. Helium density at 120 km, 300 km and 500 km versus latitude for 
"T oo =800° and the winds of Figure 29. 
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Figure 39. Helium density at 120 km, 300 km, and 500 km versus latitude 
lor Tg, = 1500° and the winds of Figure 33. 








Figure 40. Vertical and horizontal wind profiles for =800°, 50/200, 1.8. The 

values shown represent maximum amplitudes and are multiplied by sin 6 for the hori - 
zontal component and cos 9 for the vertical component, where 9 is the polar angle. 



Figure 41 . Vertical and horizontal wind profiles for = 1100°, 90/200, j3= 1.8. 
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Figure 42. Vertical and horizontal wind profiles for = 1100°, 100/200, /3 = 4.0. 
Theregion of horizontal wind labeled “negative" refers to flow from the winter hem- 
isphere toward the summer hemisphere. 
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Figure 44. Direction of wind vectors for profiles of Figure 40 
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Figure 45. Direction of wind vectors for profiles of Figure 41 
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Figure 46. Direction of wind vectors for profiles of Figure 42- 
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Figure 47. Direction of wind vectors for profiles of Figure 43. 
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that increasing the altitude of the cell (z Q ) or decreasing the wind speed at 
lower altitudes (increasing /3) generally has the effect of raising the high alti- 
tude wind speed required for a given pole-to-pole ratio. For low and moderate 
values of wind speed the variation of the logarithm of R p with W is seen to be 
linear. With increasing winds log Bp/W becomes non-linear, with the largest 
effect occurring at lower altitudes. This reflects once more the smoothing ef- 
fect of the exospheric return flow at high altitudes, and also indicates that a sig- 
nificant amount of the redistribution effect of the wind occurs in the 100 to 200 
km altitude range. 

This can be seen most clearly in the variation of helium with altitude, par- 
ticularly for an exospheric temperature of 1500° and z 0 = 180 km. Under these 
conditions, the summer pole density actually goes through a minimum with in- 
creasing altitude: The density is diminished from below by the upward wind and 
is replenished from the top by the exospheric flux. As the altitude of the 
circulation cell would most probably rise with increasing exospheric temperature, 
the occurrence of such a minimum in the density profile is not considered likely; 
however, this extreme case illustrates the result of the competition between the wind 
and the exospheric transport in influencing the vertical distribution. This low 
altitude effect is generally greater for smaller values of / 3 , i.e. when the wind 
extends to lower altitudes and when the return flow (of the wind field) is below 
the thermosphere. 

The sensitivity of the latitudinal distribution to the altitude of the wind is 
shown in Figures 48 and 49, wherea (Z Q , T^ ) (the slope of the log R p vs. w 
curve in the linear region) is plotted as a function of Z Q for the three exospheric 
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temperatures and /3 = 1.8 and 4.0. Again, the low altitude enhancement is evi- 
dent as the wind height drops, particularly at higher exospheric temperatures. 
The quantity a(Z 0 , T^) is useful also as a parameter for comparing calculated 
results with observations; given a measure of the pole-to-pole variation (e.g. 
from satellite measured densities), wind fields consistent with this variation 
can be calculated from the relation 


W= a(Z 0 , TJ log R p . 

The family of wind fields calculated in this manner for a pole-to-pole ratio of 
10 are shown in Figure 50 for /3 = 1.8 and 4.0 and the three exospheric 
temperatures. 

The variation of R p with /3 for an exospheric temperature of 1100°, W of 
100 cm/sec and z 0 of 200 km is given in Figure 51. It can be seen that for /3 
less than 2 x 10" 7 the pole-to-pole ratios at all three altitudes increase rapidly, 
with the greatest increase occurring at 120 km. This increase in R p is prin- 
cipally due to the large decrease in the helium density in the 100-150 km region 
near the summer pole, as illustrated previously for a higher exospheric temper- 
ature. The enhancement of this effect for low /3 is evident from the vertical and 
latitudinal profiles shown in Figures 52 and 53 for f3 =1.5, 2.0 and 4.0 x 10~ 7 . 
The reason for the large summer pole decrease is that upward winds in this 
region lead to an upward flux of helium which must be supported primarily by 
molecular diffusion from below the turbopause. For large winds in the lower 
thermosphere (small values of f3) this upward flux can be barely supported and 
the helium density in the 100-150 km region fells drastically. For example, for 
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Figure 50. Vertical velocity, W, required to produce R p (500 km) =10 
as function of Zg for /3 = 1.8 and 4.0. 
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Figure 52. Helium density versus alt tor T m - MU , / H 


the 100/200, T,,, = 1100°, /3 = 1.5 case illustrated, at 100 km altitude the 
helium density at the summer pole is 6.08 x 10* off*, the vertical wind is 1.7 
cm/sec, the helium scale height, H, is 32.9 km, the major gas scale height, H 1 , 
is 5.9 km and the molecular diffusion coefficient, D, is 1.32 x 10 6 cm 2 sec. 


This leads to an upward helium flux of 



= 4.7 x 10 7 / cm 2 sec. 


The maximum upward flux which can be supplied by molecular diffusion is ob- 
tained by setting the eddy diffusion coefficient equal to aero (see Johnson and 

Gottlieb, 1970); 

4> = Dn Ll + Il = 6Dn/H = 1.46 x 10 7 /cm 2 sec. 
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The difference between these fluxes must be supplied by transport down from 
the exosphere, and when this mechanism cannot provide sufficient helium the 
density falls to very low values. For the (i - 1.5 case in Figure 52 the helium 
flux is downward above 121 km, reflecting the replenishment resulting from 
exospheric transport. The sharp decrease between 90 and 120 km is due to the 
limitation on flow imposed near the turbopause: (Under these conditions the 
calculated densities become essentially meaningless. These densities result 
from differences between large numbers, where machine roundoff errors and 
numerical integration errors combine to invalidate the result). 

There is little cahnge in R p at any altitude when /3 is increased above 3. 
These high /3 wind systems are significant only at the higher altitudes 
and are characterized by a relatively strong return flow in the middle thermo- 
sphere. The variation of R p with fi shown in Figure 51 is typical of the wind 
systems studied, varying only in detail for different values of W, Z Q and T^. 

C. Comparison with Observations 
1. Satellite Data: Latitudinal Profiles 

Figure 1 gives the latitudinal distribution of helium near solstice as 
measured by the mass spectrometer flown on OGO-6 and normalized by the 
Jacchia (1965) model atmosphere to eliminate the effect of varying altitude dur- 
ing the measurement (Reber, et al., 1971). A different method of eliminating the 
altitude effect is now being employed which has the advantage of greatly reducing 
the sensitivity to the atmospheric model used. This technique utilizes only the 
exospheric temperature (from Jacchia, 1965) and the scale height corresponding 
to this temperature to extrapolate the component density to a common altitude. 
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He (crrr 3 ) at 500 km 


For helium the measurements generally lie between 400 km and 600 km altitude, 
so reducing the data to 500 km requires an extrapolation over less than half a 
scale height. Data corresponding to the same measurements as those in Figure 
1, but reduced to 500 km are shown in Figure 54; this format will be used for the 
bulk of the comparisons with the calculated distributions. 



Figure 54. Helium density measured from 0G0-6 satellite extrapolated to an altitude of 
500 km versus geographic latitude. These data correspond to those shown in Figure 1 
taken 7 June 1969 on orbit 24. 
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The solstice data referenced above indicate a density peak in the winter 
hemisphere which varies from -50° to -70° geographic latitude. This is not 
consistent with the results of the calculations presented so far, which indicate 
a cosine-like latitudinal variation, effectively mirroring the simple wind field 
assumed. Subsequent analysis of data from the same experiment (during times 
when perigee was near the poles) implies the existence of a persistent heat 
source in both polar regions, even during periods of relatively low geomagnetic 
activity (Hedin, et. al., 1970; Reber and Hedin, 1971). This postulated heat 
source is deduced from localized enhancements iii the density of molecular 
nitrogen (consistent with a temperature increase), accompanied by depletions in 
the density of helium (consistent with a rising column of air). The result of this 
polar heating is superimposed on any large scale circulation system and it 
effectively reduces the helium density in its region of influence. Thus, the direct 
comparison of the calculated helium distributions with the OGO data should be 
made with by this polar phenomenon in mind. 

The comparison of data from two OGO-6 orbits with the calculated results 
from two wind systems, chosen to match the measurements, is shown in Figure 
55. The error bars on the measurements reflect the scatter in the data, while 
the difference in location of the density peak between the two orbits is clearly 
seen. The wind cells which are characterized by different altitudes, amplitudes 
and /3’s, effectively reproduce the measurements between 70° latitude in the 
winter hemisphere to 50° in the summer. The pole-to-pole ratio at 500 km as- 
sociated with these wind systems is approximately 18; the full family of wind 
fields which yield the same R p is shown in Figure 56 for /3 = 1.8 and 4.0, and 
Z 0 between 170 and 200 km. 
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Figure 55. Data from orbits 24 and 26 of 0G0-6 extrapolated to 500 km and calculated results 
using the wind fields 270/230, /3 = 1.8 and 214/200, /3 = 4.0. An exospheric temperature of 
1100° was used in the calculation corresponding to the average daily temperature for the time 
of the measurements. 


As Z 0 increases, the distinction between the two values of /3 decreases, so 
that for Z Q greater than 230 km there is less than a 7 % difference in the high 
altitude wind speed necessary to generate the given value of R . Reference to 

P 

Figure 11 indicates that as Z Q increases the value of [i necessary to induce a 
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Figure 56. Vertical wind speed required, as a function of Z Q , to produce pole-to-pole ratio 
of 18 for helium at 500 km. This value of R p best fits the data from the OGO-6 mass spec- 
trometer. 


return flow in the thermosphere decreases, with the result that both the wind 
systems shown in Figure 55 share the common feature of a fairly intense lower 
thermospheric return flow. Conversely, lowering Z 0 and /3 while maintaining a 
given value of R p at 500 km has been shown to result in an extreme decrease in 
helium density near the summer pole in the 100-200 km altitude region, as well 
as to decrease the altitude of the return flow to below 80 km. Since rocket 
measurements do not indicate such low values for helium in the summer hemi- 
sphere, it is strongly suggested that the vertical wind profile be consistent with 
a return flow in the thermosphere. Figures 57 and 58 give the vertical profiles 
of the horizontal and vertical components of the wind systems used for the cal- 
culations shown in Figure 55. It is seen that the maximum horizontal velocity at 
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Figure 57. Vertical and horizontal wind profiles for 270/230, /3 = 1.8, T ro = 1100°. The 
region labeled negative refers to flow from the winter to the summer hemisphere. 
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the equator is about 150 m/sec for both systems; their main differences lie in 
the intensity of the return flow near 200 km and the vertical velocities below 
180 km. 

The latitudinal variation at 500 km for R p = 18, fi= 1.8 and 4.0 is shown in 
Figure 59 and 60 for a number of wind cell amplitudes and altitudes. It will be 
noticed that increasing Z 0 increases the absolute value of the helium density, 
with a 10 km change in Z Q resulting in a 10% to 25% change in density. Assum- 
ing that the eddy diffusion coefficient is as determined from equinox data, the 
absolute helium density provides a constraint on the allowable wind systems to 
explain the solstice measurements. Thus, not all the wind cells parameterized 
in Figure 56 for R p (500 km) = 18 are equally consistent with the data. 

Enhancing the effect of the vertical wind in the 100 km altitude region, 
either by decreasing Z 0 or decreasing /3 (as shown in Figure 61), results in 
lower helium densities. Winds in this region can be thought of as turbulence, 
and the reaction to a wind cell is similar to that of increasing the eddy diffusion 
coefficient: they both decrease the density of a light gas in the thermosphere. 
That this is running counter to the overall result of a wind cell whose effect is 
confined to the upper thermosphere can be seen by considering that the helium 
density at 500 km for the "no wind" case is 2.2 x 10 6 . This value is exceeded 
at the equator by as much as a factor of three for the wind cells considered 
here.. The "pumping" action taking place — transporting helium up into the 
thermosphere by the wind system - may be seen by reference to Figures 62 
and 63, the latitudinal and vertical profiles of helium associated with a circula- 
tion cell when only the amplitude of the wind is varied. Increasing the wind speed 
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Figure 59. Helium density at 500 km versus latitude for = 1 100°, /3 = 1 .8 
and 215/210, 270/230, and 365/260. These wind systems all produce nearly 
the same R p (500 km), but the absolute values differ by more than a factor of 
two. 
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Figure 60. Helium density at 500 km versus latitude for = 1100° = 4.0 

and 176/180, 214/200, 280/230, and 380/260. 
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Figure 61. Helium density at 500 km versus latitude for = 1100°, 
260/230, and y8 = 1.5, 1.7, and 4.0. 
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Figure 62. Helium density at 500 km versus latitude for = 1100°, 
/S = 4.0, and 200/230, 260/230, and 300/230. 
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Figure 63. Helium density as a function of altitude for the same conditions as Figure 62. 

primarily increases the density in the winter hemisphere while the density 
in the summer hemisphere is relatively unchanged. The diffusion limitation 
exhibits itself in both hemispheres: the vertical profiles at the summer 
pole are nearly identical up to 150 km, while at the winter pole there is a 
"piling up" of helium which is transported in by the wind and cannot readily 
diffuse down below 150 km. 

1. Rocket Data: Vertical Profiles 

One of the two main conclusions from the many determinations of helium 
density by rocket-borne mass spectrometers is that the lower thermosphere 
exhibits a similar seasonal variation to that observed at higher altitudes from 
satellites. The other, perhaps more significant, result is that helium often is 
not in static diffusive equilibrium with file major gases in the altitude region 
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from 120 to 250 km. On the contrary, in nearly all the measurements reported 
to date by the group at the University of Minnesota (e.g. Hedin and Nier, 1966; 
Krankowski, et al., 1968) and one by Goddard Space Flight Center (Cooley and 
Reber, 1969) the altitude profile of helium has indicated a lower scale height, 
Hjjg.than would be consistent with the temperature deduced from the scale heights 
of the other gases. The single exception is the winter measurement at Fort 
Churchill reported by Hartmann, et. al. (1968; see also Muller and Hartmann, 
1969) which indicated a high density and a nearly static scale height. A number 
of the Minnesota results have been summarized and interpreted by Kasprzak 
(1969) as due to an upward flux of helium, perhaps initiated by lateral transport 
in the exosphere (McAfee, 1967). Reber (1968) suggested that the relatively long 
diffusion times in the lower thermosphere coupled with a temperature change or 
variation in turbopause level might be responsible for the low scale heights. 

None of the mechanisms proposed, however, indicate why the flux is predomi- 
nantly upwards, or equivalently, why the scale heights are generally lower than 
expected. 

Reference to Figures 64 through 68 indicates the behavior of the scale 
height of helium under the influence of a variety of wind cells. Figure 64 shows 
the summer profiles up to 500 km for fi from 1.5 to 4.0 compared with the static 
profile; it is seen that all the scale heights approach the static value at high 
altitude, even though there is as much as a factor of two difference between 200 
and 300 km. Figures 65 and 66 emphasize the lower thermosphere summer and 
winter profiles for the same family of wind systems, while Figures 67 and 68 
give similar profiles for a set of wind systems (both for /3 = 4.0) which produces 
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Figure 64. Helium scale height, H^ e , as a function of altitude at the summer pole for 
To, = 1100°. The winds represented are 260/230, = 1.5, 1.7 and 4.0 ; al so shown i s 

the scale height in the case of no wind. 
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Figure 65. Same as Figure 62 with emphasis on the region below 300 km. 


pole-to-pole latitudinal ratios consistent with satellite measurements, but which 
have different altitudes and amplitudes. 

The dynamic summer profiles exhibit a lower scale height than the static 
scale height for /3 1 1.7 and altitudes less than 170 km. As the value of /3 is 
increased (raising the effective altitude of the cell), the altitude of lower scale 
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height is also raised. For lower altitude cells (and for altitudes above about 


200 km) the downward flux from the exosphere dominates the distribution and 


the dynamic scale heights are greater than their static counterpart. For /3 = 


1.5 the dynamic scale height is less than the static up to 120 km, indicating the 


region of dominance for this particular cell. 
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Figure 67. Helium summer pole scale heights for = 1100°, /3 = 4.0, 176/180 and 380/260, 
emphasizing the result of lowering the dominant altitude of the wind field. 

The wintertime profiles show that essentially the opposite situation pre- 
vails in a region of subsidence. Here, the cells effective to lower altitudes 
produce scale heights lower than static in the lower thermosphere, while higher 
altitude cells generate higher scale heights. At the upper altitudes all the 
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300 



dynamic scale heights fall short of the static values due, once more, to the up- 
ward flux to the exosphere. This behavior is shown accentuated in Figure 69 
where the scale height due to a very low altitude cell (/3 = 1.8, Z Q =170 km) is 
compared to those from higher altitude cells and to the static scale height. In 
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Figure 69. Helium winter pole scale heights for = 1 100°, fH = 1.8, 
365/260, 215/210, and 125/170. 


this case H He is nearly a factor of two lower than static in the altitude region 
from which most of the rocket data are obtained. Thus, while virtually any of 
the winter wind profiles considered generate helium profiles consistent with 
rocket measurements, these same measurements require a summer wind systexr 
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which is more effective in the lower thermosphere. It is most likely that the 
simple, symmetric wind system considered here is not adequate, and also that 
there are other cells (e.g. diurnal) existing in conjunction with the seasonal cells, 
each independently influencing the helium distribution. 

D. Time Development of Response 

The evolution of the helium response to a wind system that is suddenly 
"turned on" was studied using time increments of three hours, twelve hours and 
two days. This was found necessary as use of large increments caused oscil- 
lations in the early phases of the response, while use of small increments re- 
quired excessive computer time. The combined results are shown in Figures 
70 through 73 for exospheric temperatures of 800°, 1100°, and 1500° using /3 =1.8, 
and 1100° for (3 =4.0. Here the helium densities at the winter and summer 
poles are given as functions of time for altitudes of 120, 300 and 500 km. It 
can be seen that there is an initial response which varies from about fifteen 
days for an 800° exosphere down to three days at 1500°, followed by a relatively 
long term density variation. In the case of all three temperatures, the latter 
manifests itself as a gentle increase (less than 0.6%/day) at both poles and all 
altitudes. It is also apparent that the wind is more effective in evoking a re- 
sponse at the higher altitudes as the initial phase is significantly longer at 120 
km for all three temperatures, indicating that the variation is being propagated 
downwards. 

The exact cause of the long term portion of the response is not clear, but it 
is apparently related to the exospheric transport (through the induced "pumping" 
effect discussed previously) and the relatively long time required for helium to 
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ure 70. Time development of the summer and winter pole helium distributions at 120, 300 
and 500 km for low solar conditions; the 70/120, /3 = 1.8 wind is “turned on” at t = 0. 


se up through the lower thermosphere. This phase of the response is 
ically less important than the initial phase, however, as it would most 
y be masked by shorter term variations in exospheric temperatures due to 
letic activity and the 27-day solar cycle. 
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Figure 72- Time development of the helium response for medium solar conditions 

and a 214/200, (3 = 4.0 wind field. 


enhancement of a factor of four over the density at the sub-solar point. The 
amplitude drops below this for both high and low solar cycle conditions; however, 
the temperature differential at high solar cycle would tend to produce a similar 
variation due to exospheric transport (McAfee, 1965), so the net effect might be 
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equal to or larger than that at mid-solar 
pears likely that less than a factor of two 
at high altitudes due to this mechanism. 


cycle. At low solar conditions, it ap- 
day-night variation could be maintained 
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E. Other Minor Species: Argon 

Of the minor gases of interest after helium, argon is the most useful to study 
as it is also inert and it has been measured by rocket and satellite borne mass 
spectrometers. Also, its mass of 40 is greater than the mean mass in the 
lower thermosphere, so in accord with Equation 5 its response to a wind system 
should be opposite that of helium. In addition, due to its high mass (relative 
to helium) the effect of exospheric transport should be negligible. 

The latitudinal variation of argon at 300 km (where it can, in principle, be 
measured by satellite-borne mass spectrometers) is shown in Figure 74 for the 
214/200, /S = 4.0 and 270/230, = 1.8 wind systems which provided good agree- 

ment with helium measurements at 500 km (see Figure 55). The effect of the 
high relative mass can be seen immediately, with the densities near the summer 
pole higher than the winter densities by nearly a factor of four; the distributions 
corresponding to the two systems also compare well with each other, with less 
than a 4% density difference at any latitude. As the effective altitude of the 
vertical wind is raised, by raising /3 (Figure 75) or by increasing Z 0 (Figure 
76), the effect on argon is to decrease the amplitude of the latitudinal variation 
in a nearly symmetric manner. It will be recalled that the same variation in 
wind field caused a general increase of the helium density at 500 km, while 
maintaining the pole-to-pole ratio relatively constant (see Figures 60 and 61). 
Thus, when simultaneous latitudinal profiles of argon and helium are available, 
one can in principle narrow down the family of wind fields consistent with the 
two distributions. 
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Figure 74. Argon density at 300 km versus latitude for T ro = 1100°, 214/200, /3 = 4.0 and 
270/230, /3 = 1.8 winds. These winds give the best fit to the OGO-6 data for helium. 
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Figure 75. Argon density at 300 km versus latitude for T 
260/230 and fl = 1.5, 1.7 and 4.0. 
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The insensitivity to exospheric transport can be seen in Figure 77 which 
shows vertical profiles for an exospheric temperature of 1500° and vertical 
wind speeds of 4m/see. Under the same conditions (see Figure 33), helium 
displays a distinct depression in the summer hemisphere near 200 km, with 
the higher altitude density enhanced by exospheric flow. 

IV. CONCLUSIONS 

A large scale meridional circulation system in the thermosphere (upwelling 
in the summer hemisphere, flowing toward and descending in the winter hemi- 
sphere) was shown to be sufficient to generate the observed enhancement of 
helium in the winter upper atmosphere. The increase of exospheric transport 
with temperature results in a smaller latitudinal variation at high solar activity 
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than at low activity, due to the large smoothing effect of the return flow. 
Horizontal diffusion in the thermosphere, however, is negligible as a smoothing 
agent compared to exospheric flow. On the basis of satellite-type measure- 
ments of helium alone it is impossible to distinguish between a variety of wind 
fields as the causative mechanism; however, wind fields consistent with the 
helium distribution measured by OGO-6 are characterized by vertical velocities 
of two to three meters per second above 200 km and horizontal velocities at 
the equator of one to two hundred meters per second. These are within a factor 
of two of the amplitudes proposed by Johnson and Gottlieb (1970) to explain the 
temperature in the winter thermosphere, but are 100 km higher in altitude. 

Argon is affected in the opposite way from helium, being enhanced in the sum 
mer hemisphere and depleted in the winter; there is negigible effect here from 
exospheric transport. The calculated vertical helium profiles indicate departure 
from a static-diffusion profile in much the same manner as observed by rocket 
measurements. In order to be more consistent with observations, however, it 
would be necessary to disregard the simple, symmetric circulation cells used 
here and adopt a wind field which is effective to lower altitudes in the winter hemi- 
sphere than in the summer. By use of latitudinal data on helium and argon, in 
conjunction with vertical profiles of helium, it should prove possible to narrow 
down the number of potential wind fields causing the distributions. Disregarding 
photochemical effects, atomic oxygen should exhibit the same behavior as helium, 
but withalower amplitude. However, since it is more temperature sensitive than 
helium, the wind and temperature effects tend to be self cancelling. Thus, the 
net effect in oxygen might well be the absence of an expected enhancement in 
the summer hemisphere at higher altitudes. This type of response has already 
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been noted in the high latitude neutral gas data from OGO-6 following magnetic 
storms (Taeusch, et. al., 1971), where the N 2 density rises, the helium falls 
and the oxygen remains relatively constant. 

The time response of the helium density to a wind field indicates a signifi- 
cant variation in less than half a day. This leads to the likelihood of a factor of 
two to four density enhancement at night. This, as well as the other effects 
discussed here, clearly indicate the need for the inclusion of dynamics in 
describing and studying any but the simplest of upper atmosphere phenomena. 
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APPENDIX A 

COUPLED MOMENTUM AND CONTINUITY EQUATIONS FOR A MINOR GAS 
The continuity equation for a single gas is written as 


|_^.+ V- (n v) = 0 


where n is the gas density in molecules (or atoms) per cubic centimeter and v 
is the flow velocity of the gas. In spherical coordinates, with no longitudinal 
dependence, this becomes 


T7 <" + 


2 n v 


1 B , , m 3n 

— - ( n v„ sin 6 )= , 

r sin 9 B 0 6 ’ B t 


(A.l) 


where 0 and r are the polar angle (latitudinal) and radial variables, and the 
subscripts refer to the respective velocity components. The momentum equation 
for a neutral atmospheric component experiencing negligible acceleration, can 
be written 


Vp_ -nmg+mnv 

n 


[v - V] • = 0, 


where 

p n = partial pressure of species with density n, 
g = local acceleration of gravity, 
m = molecular mass (gms) of species n, 

v = momentum transfer collision frequency for gas n in background gas, 
and 

V = flow velocity of background gas. 
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Using the ideal gas law, p = nkT, this becomes 


or 


where 


n [v - V] = [n m g - Vp ] 

m v n 


n [v - V] 


kTn 


m v 



m g 
k T 


k = Boltzmann constant 
and T = local temperature. 

The radial component of the momentum equation then becomes 


n [v r - V r ] = - D 


3 n n (1 + a) 3 T n 
+ — 1 ' + — 

3 r T 3 r H 


(A.2) 


where 


H = 


D = 


kT 

mg 

kT 

mv 


, the local scale height for the species with mass m, 
, the molecular diffusion coefficient, 


and a is the thermal diffusion factor (Chapman and Cowling, 1939; Kockarts and 
Nicolet, 1963; Colegrove, et. al., 1966). Similarly, the latitudinal component 
becomes 


n[v e -VJ =-2 [I" "»+<>) II 
99 r [_3 0 T 30 


(A .3) 


Lettau (1951) has rigorously modified the atmospheric diffusion equation to 
include the effects of eddy diffusion. Colegrove, et. al. (1965) arrive at the 
same expression for the vertical concentration gradient as Lettau by consider- 
ing the flux for a given component to be composed of a molecular diffusion term 
and an eddy diffusion term. Following this approach, the radial momentum 
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equation becomes 


n [v f - V r ] = - D 


3 n n(l+a) 3 T n 
3r + T 3 r + H 


- K 


Bn n 3 T n 
BT + T T7 + H 7 


(A .2 ') 


where K is the eddy diffusion coefficient and H' is the scale height associated 
with the mean molecular mass. No corresponding eddy mixing term is added 
to the latitudinal flux equation as no horizontal temperature gradient is assumed, 
and the horizontal density gradient for the minor gas is assumed to be small in 
the region where eddy mixing is appreciable. (Since the addition of this term 
would decrease any horizontal gradients, the results of the calculation when it 
is neglected justify the latter assumption.) 

Adding -B/Br (n V r ) to both sides of (A.l), the continuity equation becomes 


_n(v r -V t ) = - 


2 n v 


- — (n v fl sin - — (nv ) 

r sin 9 B 9 K 9 ' Br 


B n 

b7' 


(A.l 1 ) 


adding 


1 

r sin 9 


Yq ( n V* Sin 9) 


to both sides of (A.l’) and making use of (A.3): 


B 

B r 


n (v r - V r ) - 


1 

r sin 9 


B 

b e 


(n V e s in 9) - 


2 n v f 
r 


BV r 

n 

B r 


- V. 


B n 

B r 


+ 1 JL flsin e u) 

r sin 0 30 [_r ^3(9 T B 9 J 


(A .4) 

B n 

bT' 


The continuity equation for the major background gas with number density N is 
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3 2 N V r I 

W /»T TT \ r J. 


3 r 


(N V r ) + 


r sin 8 3 6 


3 (N Y, sin d) = Ail 
; B t 


Assuming no change with time or latitude for the major species 




this can be written 


r sin 6 B 9 


<y 6 


1 -S 2 V 

sin 0) = -A —— (N V r ) + r 

N 3 r r r 


9 3 N (A.5) 

N r B 9 ' 


Using the identity 


B n 

B 6 


substitution of (A.5) into (A .4) yields 


T-r " <v - - v .) + ¥ r- <N v ,) + — * — — -- — - ^ 

3 r NBr r Nr B0 r B 6* r 


(A .6) 


3 V 

- n -_L - V _ + 

r3r r 2 sin0 30 


■r ..BN 1 


n sin 8 (A! 1 + n (l +<X) Bt \ _ 3 n 
V 6 T B e)\ B t 


Replacing n (v f - V f ) by use of (A.2') and rearr anging gives 

A? - 3 Jn T 3 n , n Cl + a) 3 T , n"| v [b n , n 3 T . nil 

ifTFrljT* T — 3T hJ [ar + T r? + H]J 

♦ i (d b U + z\ + K [A + fi il + J] 

r l L 3 r T B r hJ L 3 r T B r H'JJ 

+ v [fL AA - 3 n l , y ^ F n 3 ^ 3 nl 

r (n B r B r J 9 r j^N 3 9 ~ B ej 

+ I AL [psingfAil ^Cl + a) Al\l. 

r 2 sin^ 3 ^L \ 3 ^ T 3 ^/J 


(A. 7) 
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This is the continuity equation for a minor gas, modified by motion in the 
background gas, which appears as Equation (4) in Chapter II. 
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APPENDIX B 

f 

MODEL ATMOSPHERE 

The model atmosphere used in the calculations was based in the COSPAR 
International Reference Atmosphere (CIRA, 1965) for the altitude range 80 to 
120 kilometers and the Jacchia (1965) model as modified by Walker (1965) for 
altitudes above 120 km. The CIRA model is presented as a tabulation and util- 
izes a number of straight line temperature profiles. In the calculation of hori- 
zontal winds, the expression for B (r) contains a term proportional to the scale 
height of the major species, which in turn, is related to the temperature 
gradient: 


3 N N B T N A 
Tf + 7 Tr + H 7 = °‘ 


It was found that the horizontal wind so calculated went through a number 
of discontinuities at the intersections of the straight line temperature profiles, 
so the tabulated temperature profile was modified slightly to eliminate the dis- 
continuities in slope. The CIRA and the modified temperature profiles are 
shown in Figure 79; the effect on the calculation of B (r) for both profiles is 
shown in Figure 80 for a typical wind system. The smoothed temperature is 
given in Table B1 and the complete tabulation, including the densities and 
mean mass up to 120 km, is given in one of the block data subroutines listed 
in Appendix E. 

Above 120 km the model is analytic and presents the temperature, T, and 
component number densities, n., as functions of altitude, z: 
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Figure 79. CIRA, 1965 temperature profile compared with the smoothed 
profile used for the present work. 


T ( z ) = ^ - (To, - T ) exp (- a <f ) 


(B.l) 


and 


n i (z) = n. (120) 


1 - a 


1 - a exp (- ct £) 


1 +a+y 


exp(-a-yg), (B.2) 


where 

= exospheric temperature, 

T 12o = temperature at 120 km, 
a = S + 0.00015, 

S = 0.0291 exp (_ x 2 /2), 

€ = (2 - 120) (R + 120)/R + Z = geopotential altitude, 
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Figure 80. Effect on (twice horizontal wind component) of smoothing 
CIRA 1965 temperature profile. 
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Table B-l 


Altitude 

(km) 

CIRA 

Smoothed 

CIRA 

Altitude 

(km) 

CIRA 

Smoothed 

CIRA 

80 

186.0 


101 

212.2 

212.6 

81 

186.0 

184.6 

102 

215.7 

216.0 

82 

186.0 

184.6 

103 

220.0 

219.6 

83 

186.0 

184.7 

104 

224.6 

223.4 

84 

185.9 

185.0 

105 

229.0 

227.6 

85 

185.9 

185.3 

106 

233.4 

232.1 

86 

185.9 

185.8 

107 

237.9 

236.9 

87 

185.9 

186.5 

108 

242.3 

242.0 

88 

185.9 

187.3 

109 

246.8 

247.6 

89 

185.8 

188.3 

110 

251.1 

253.5 

90 

185.8 

189.3 

111 

261.6 

260.0 

91 

188.4 

190.6 

112 

271.9 

267.0 

92 

190.9 

192.0 

113 

282.3 

274.5 

93 

193.5 

193.6 

114 

292.7 

282.6 

94 

195.9 

195.3 

115 

302.9 

291.5 

95 

198.2 

197.2 

116 

313.1 

301.3 

96 

200.4 

199.2 

117 

323.6 

312.2 

97 

202.4 

201.5 

118 

334.0 

324.5 

98 

204.4 

204.0 

119 

344.0 

339.1 

99 

206.3 

206.7 

120 

355.0 

355.0 

100 

208.1 

209.5 
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X = - 800/750 + 1.722 x 10' 4 (T^ - 800) 2 

R = radius of earth = 6356.77 km, 

a — f(n — ^'l20^'^' 00 

a = thermal diffusion factor (0 for all but helium) 

y = m i Bi2o /cr k T.. 

g i 2 o = acce l era tion of gravity at 120 km = 944.655 cm/sec 2 
and 

k = Boltzmann's constant. 

At 120 km the temperature is 355 °K and the number densities are: 

n (N 2 ) = 4.0 x 10 11 cm -3 
n (0 2 ) = 7.5 x 10 10 

n (0) = 7.6 x 10 10 . 
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APPENDIX C PRECEDING PAGE BLANK NOT FILMED 
RELATIONSHIP OF HORIZONTAL AND VERTICAL WINDS 

The horizontal component of the wind field is related to the vertical com- 
ponent through the continuity equation for the major species (assuming 3 N/3t 
= 0 ): 


3 

3 r 


(N V r ) 


2 N V r 

+ 

r 


+ t — t A (N V 5 sin 9) = 0, 

r sin odd ° 


(C.l) 


where 

N = major gas number density (sum of O, 0 2 , N 2 ), 
V r = radial component of wind field, and 
\ g = latitudinal component of wind field. 
Rearranging and noting that 


1 3 N _ J_ 
N 3T ” H 7 ’ 


3 

3 6 


(N W 9 sin 9) = - t sin 6 — — (N V ) - 2 N V sin 9 

dr r 


= - N sin 9 




+ 2 V 



(C.l*) 


Using the expansion of the vertical wind component, 


V r ('• 0 = 21 \ < r > P ^ (")• 

t 

and the assumption of no latitudinal variation in N, (C.l') becomes 
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where 


d 6 


(sin 6 \ g 


(r)) = - sin 6 



( r ) Pfc (M). 


(o = 


/ 3 \ (O 

l r -T7"- 


r Vo 
H' (r) 



With fj.- cos Q and d/j. = -sin ddd, we obtain 


Integrating over *± from /j.’ to 1 (5' to 0) leads to 


V, = -£ - fpV>d M 

V <w 2 ) 172 J»» * 


P^ 1 (^) (Magnus and Oberhet inger, 1949) (C.2) 

l 


where 


? ? (M) = 1 f Pp (M) d M 

(1-V 1/2 * 


r(0 p j (m) 
rV + 2) 


(C.3) 


Vi 00 - u Pj (fx) 

(t + 1 ) (1 - M 2 ) 172 
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APPENDIX D 


METHOD OF SOLUTION 


1. Harmonic expansion 

The minor gas continuity equation modified to include motion in the back- 
ground gas was shown in Appendix A to be : 


3 n 

Ft 




4 “ 


3 n n ( 1 + a) 3 T ^ r» 
3~7 f T 3~7 H 


+ K 


3 n n 3 T n 
Tr + T 3~r + H 7 " 


+ V. 


n 3 N 3 n 

N 3 r 3 r 


+ V *7 


n 3 N 3 n 
N 34? “ Ye 


1 


r J sin 0 ^ L 


a /3 n n (1 + a) 3 T' 

D s in 6 — - + — t { ___ 

V3 <9 T 3 e t 


(D.l) 


It was assumed in the solution of (D.l) that neither the temperature, T, nor the 
major component number densities, N, varied with latitude; this implies also 
that H, H', and D are 8 independent. The minor gas number density n(r,£) was 
expanded in a series of Legendre polynomials 

n (r, 8) = 2 ^ n n ( r > P n (F>- (D.2) 

n 

where 

ijl = cos e, 

and a solution was sought for the n th coefficient, n n (r). The horizontal and 
vertical components of the wind were also expanded in Legendre series and the 
full wind field was expressed in terms of the coefficients for the vertical com- 
ponent (see Appendix C): 
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and 


V r ( r - *) = 2] 

l 

v e ( r - d ) = - Y1 b * (m) 
* 


(0 = 


3 r 


+ r — 


N 


3 N 
B r 


+ 2 \‘ 


After the above simplifications and substitutions each term in (D.l) is 
multiplied by (/j.) and integrated from -1 to +1 (9 = 0 to 6 = tt, as a negative 
sign has come in through d/x = - sin 9d 9). Thus, the equation for the coefficient 
of the m th harmonic becomes 


2 _ B f n / 3n m n n ,( 1 + a ) 3T ( 9 "n, ^ 3 T 

2m+lBt 3 r ^ \ 3 r T Br H / \ 3r TBr 


. Oil 2 
+ H'/J 2 m + 1 



+ 


n m C 1 + a ) 

T 


_3 T 
B 


T n m \ „ 
- + — ) + K 

r H / 



n x n V| 2 8 

mol mil nm 

+ TBT + H 7 /J 2m+ 1 



\ 


n Bn 

n n 




(M) P m (M) d M 


(D.3) 




+ 1 , Bp (fj.) 

p^ 1 (M) P„ (m) dM 


B Q 


BY -1 f 1 3 / . B p n 0)\ 

ML” 1 sin~9 m?{ Sln$ ^uT-j =°- 

n ^ 


Here 


1 = J_ 3_T 1 _ 1 BN 

TBr + H' NBr 
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With the substitutions 


and 


Atnm = f P n M P m (/*) d M- 

•'-1 

r +i bp (m> 

% nm = J -TT~ P * (M) d * 



B 

B 0 


s in <9 


3 P n fa) 

b e 


d/u. 


equation (D.3) is equivalent to (10) in section EC: 


B n 


2 m + 1 Bt 2 m + 1 B 


2 m + 


i li }-ZZ v -t 


•C,„ 


n Bn 

n n 

+ 


~T L Vn B L 4 ^ n " 


where 




c , 

nm 


H* 


B r 


(D.3*) 


= D 


' B n r„ n m ( ■*■ + a ) B T n m' 

— + 1 + 

B r T B r H 


+ K 


B n n rp 
m mol 

TT + T ¥7 + 


H*J 


2. Numerical integration 

a. Lindzen and Kuo algorithm 

A numerical solution to (D.3 1 ) was obtained by use of an integration technique 
described by Lindzen and Kuo (1969). They express a differential equation of 
the form 


d 2 f 
d x 2 


, N d f . 
+ g (x)— - + h 
d x 


(X) f = r (x) 


as the finite difference equation 


(D.4) 
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(D.5) 


where 


A. f iM +b 4 f 4 +q f l+1 =D if 


A 1 g (Xi) 

1 - v2 2 S x 


(Sx)= 


B, = - ■— - + h(x. ), 
(S x)2 


c. =_L_ +1^ 

1 - v2 2Sx 


(S x y 


Di = r ( Xi ), 

Sx is the finite-difference grid interval and i = 1, 2, 3, 
conditions 


d f 


and 


become 


and 


+ a f = b, at x = 0 

d x 1 1 


d f 


, + a, f = b, at x=l 

dx : * 


K f o + % f, = D„ 
A, f,., +B, f, =D,. 


The difference equation is solved by substituting 


f. , = a. . f. + ( 3 . , 

1” 1 1-1 1 'l-l 


(D.6a) 


(D.6b) 


(D.6c) 


(D.6d) 


.1-1. The boundary 


(D.7) 
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into it and obtaining 


a. = . 


-C r 

A i a i-l + B i 


and 


Pi = 


Dj - A i A-i 

A i a i-i + B i 


The lower boundary condition becomes 


B, D. 

a = - — and /3 = 

° \ ° A b 


Thus, knowledge of f f provides all the f . through (D.7); f x may be found by 
substituting (D.7) into the top boundary: 


f i = 


D t ~ A t Pj-i 

B t + a i-i A t 


b. Time dependent solution 

The one dimensional solution to (D.3’) is obtained by expressing it in the form 
of (D.5)/by use of the finite difference approximations (Crank and Nicolson, 1947) 

B 2 n^ n i + i - 2n i + n i-i 
B r 2 (S r) 2 

B n n i + i ~ n i-i 
B r 2 8 x 

3n H n | + 1 - "{ 

B t St 

With these substitutions it becomes 
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d +1 - n} a. 
1 11 


A t 


n l+i - 2 n] + nj_i nj;; - 2 n{*‘ + nj!' 
+ 


(8 r) 2 

(8 r) 

"l-i n ltl - 

c. 

1 r 


8 r 


8 r 


where a if b A and c £ are the coefficients in (D.3*). Rearranging: 


n j + i 
i+i 


r a. b. n 

i i 

+ nj + 1 

3 

2 a. "I 

1 4- c. 

+ rd + 1 

"a. b. 

i i 

_(S r) 2 ' 2 8 r . 

8 t 

(8 r) 2 l J 

1“1 

_(S r) 2 2 8 r _ 


(D.8) 


' a i b i ' 

— n) 

2 2 a i 

— 4- r 

— 

a. b. ' 

(8 r) 2 2 8 r _ 

i 

* (S o 2 ‘J 


L(S r ) 2 2 8 r J 


Consistency with (D.4) requires a coefficient of unity for the second derivative so 
we divide by a i . Then, putting the result in the form of (D.5) yields 


. _ i V a i 

Aj , 

(8 r) 2 2 5 r 


B. = 


2 c i 2 




and 


1 (8 r) 2 a i a i 5 t 

r 1 V a i 

i_ (Sr) 2 + 2 Sr’ 

D 1 =-A 1 n{. 1 -^ i+ _l T ^n{-C 1 n{„ 


(D.9a) 


(D.9b) 


(D.9c) 


(D.9d) 


for the coefficients . 

The general solution to (D.3 1 ) is then obtained by straightforward extension 
of this technique utilizing an L-dimensional vector as dependent variable (L 
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corresponds to the number of Legendre polynomials used in the expansion) and 


matrix coefficients. (D.3*) is rewritten 


where 


L 

L 


N=1 L. 


B 2 B B 

°M , N ~ ~ + ^M,N 3 r + C M,N ^ 


B 


= 0 


®M,N ^MN 


M,N 


3 (D) + D . ± K 

V ' tl IT' 


3 r H H' 


o 2 M + 1 \ 1 .1 1 

MN ” 2 4—1 ^MN[ D + K 


“M,N 


B /d\ „ B / l 1 

Br H +K Br VH'J 


s 2M+1 DC iin 

§MN+ — — 


2 M + 1 Y' \ V l . i I I i 

2 Z_i Lh' + 7 \ ^ 

1 = l 


+ K 


Then the coefficients, (D.9) become 


A s 


MN 


MN 


(S r )2 2 S r 


B. 


2 8 , 


MN 


(«r)5 


C . ... 


2 S, 


MN 


MN 


(D + K) S t 


C 


£ L 

MN MN 

+ 


(S r) 


2 2 8 r 


^MN 

^Mn 

nl t - 

^ 2S mn 

p 2 ^MN 

4- 1 . 4- 

(S r) 2 

28 r 

i-l 

_(S r) 2 

* N (D + K)8t 


S MN ^MN 

(S r) 2 2 S r 


rv 


i + 1 


(D.lOa) 


(D.lOb) 


(D.lOc) 


(D.lla) 


(D.llb) 


(D.llc) 


(D.lld) 
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The steady state solution can be obtained by dropping the third term in B t and 
setting Dj equal to zero. 

c. Evaluation of A 9 , B p , and C 

-1MN O/MN MN 

The coefficient A .0 is calculated directly 

'V mn 

A -C«N = j F i <U) P » (U) P " ^ d U 


in the course of the machine integration of (D.3'), using a machine supplied 
subroutine to perform the integration over the appropriate interval and program 
supplied polynomials. 

The second coefficient, 

r +1 3P n (u) 

B -e MN = J ^ -fe- p « d u 

•'-1 

is reduced to a tractable form by the substitutions 


3 P„ (u) 

a e 


s in 6 


* F n ( U > 

3 u 


= - (1 - u 2 ) 1 ' 2 


3 u 


J_ p (u) _ " < UP n ( U > - P n -1 (»)) 
BU " _ ( U 2 - 1) 


r t' M = -iW) P i w 


= -<l -u 2 )" 2 JLp{(u> 
(u (u) - P t _ i (u)) 

(i + 1) (1 - U 2 ) 1/2 
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Thus, reduces to 



n [P e _ i (u) - u (u)] • [u P n (u) - P n _ x (u)] P ra (u) d u 
(l + 1) (u 2 - 1) 


This integration is carried out by machine, with the singularities at ± 1.0 being 
avoided by using limits of ±0.99. 

The coefficient 


*-mn 


p P-W a 

J sin 6 3 6 


is simplified by the observation that 


1 J • P n < U > 
sin 6 


sin 8 3 6 


3 6 


L sin e -^0 — J du 

j + n (n + 1) P n (u) = 0 


from Legendre's equation. Thus 


= _ n (n + 1) 



(u) P n (u) 


d u 


2 m (m + 1) 
2 m + 1 


for n = m, 


- 0 for n/m. 
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PRECEDING PAGE BLANK NOT FILMED 

APPENDIX E 

PROGRAM USED FOR SOLUTION OF 
MINOR GAS CONTINUITY EQUATION 


IMPLICIT REAL*8(A-H,0-Z) 

COMMON/INDEX/L, M, N, I, LI* Ml, N1 
COM MON/CALC/DR* CS, RE, ALTO, TINF, EDC , TLB, ZLB, MN2, M02, 

1 MOl , MHE 

COMMON/MOD/TEMP ( 41 ) , MM(41), DENN2<41), DEN02<41), DEN0K41), 

1 DENHE ( 41 ) , DEiMA ( 41 ) 

R E A L * 8 

1 MASS, MM, MN2 , M02 , MOl, MHE, MBAR, NUM 

COMMON/COEF/ ALMN ( 6 , 6 , 6 ) ,' BLMN(6,6,6) 

DIMENSION XD EN ( 6 , 6 ) , DENSI422), ZET1I422), ZET21422), ZET3I422) 
DIMENSION UDEN (6,6), BL ( 6 ) , BP ( 6 , 6 ! , DM 1 ( 6 ), , DM2 ( 6 ) , DM3 ( 6 ) , DM4 ( 6 ) 
DIMENSION LL ( 6 ) , MV(6) 

DIMENSION AL ( 6,6, 422 ) , BE(6,422), F(6,421), VL(6), A(6,6 ), 

1 B ( 6 , 6 ), C ( 6 , 6 ), DEN (6,6 ), PR0(6,6 ),PR0B(6 ), DM ( 6 ) , 

2 NUM ( 6 ), ALF ( 6 ), HEIGHT(422) ,FM(6) , AL P ( 6 , 6 ) , PN ( 6 ) , DPN ( 6 ) 
DIMENSION AUX ( 200 ) ,DH1( 422 ) ,DH2( 422 ) ,DH3( 422 ) 

DIMENSION SCHTK422), SCHTA(422), DIF(422), SHI1422), SC I ( 422 ) 
DIMENSION VLK 422), BLII422), DDHI ( 422 ) , EDC I ( 422 ) 

EXTERNAL PLP,PLL 

DATA AL/9 *0 . / , BE / 1 . 663 E 09 , 2*0 . / , TD/-0 .4/ 

DENT ( ALT ) = DN2 ( ALT ) + DOKALT) + D02 ( ALT ) 

SCHT ( ALT , MASS ) = 8 . 3 1 E07*T ( ALT ) * ( < RE + ALT ) **2 > / l MASS*980 . 665* ( RE 
1 ***2)) 

DI FC( ALT )=( 1 .69E1 9/DENT ( ALT ) )*< ( T ( ALT ) / 27 3 . 16 ) **0 . 69 1 ) 

Y ( ALT) = 1 . /( D I FC ( ALT ) + EDC) 

DT(ALT)=(T(ALT + DR) - T ( ALT-DR ))/( 2 . *DR ) 

TDT< ALT)=DT(ALT)/T( ALT) 

SH ( ALT ) = ( 1 . + TD)-TOT(ALT)+ 1 ./ SCHT ( ALT , MH E ) 

SC ( ALT ) =TDT ( ALT ) + 1 . / SCHT ( AL T , MB AR ( ALT ) ) 

DD( ALT)=(DIFC( ALT + DR) - DIFCIALT - DR)) /(2.*0R) 

DDH(ALT) = (DIFCIALT + DR) * SH (ALT + DR) - DIFC 

1 (ALT - DR) * SH (ALT -DR) ) /(2.*DR> 

DHH(ALT) = (SC (ALT + DR) - SC 

1 (ALT - DR) ) / ( 2 . *DR ) 

DHE ( TH , I ) = F ( 1 , I ) + F < 2 , I )*DC0S(TH) + 0.5*F(3,I)* 

1 (3.=:= DCOS(TH) * DCOS(TH) -1.) 

2 + 0.5* ( 5.*DC0S( TH)** 3-3 . *DCOS < TH > ) *F < 4, I ) ' 

3 + (1./8.)* < 35.*DC0S(TH)**4-30.*DC0S(TH)**2‘ +3. >*F( 5, I ) 

4 + (1./8.)* (63.* DCOS(TH) **5-70. *DCOS(TH)* *3+15. *DCOS(TH)) 

5 * F ( 6 , I ) 

5 FORMAT ( 1 1 • ) 

500 FORMAT! 15, 1P7E16.6) 

250 FORMAT ( 3X, 1P6E16.6 ) 

251 FORMAT ( T 1 20 , ' A 1 ) 

252 FORMAT ( T 1 20 , 1 B 1 ) 

253 FORMAT LT120 » ' C ' ) 

50 FORMAT! 3X, 1P6E16.6 ) 

51 FORMAT ( T 1 20 , 'DEN') 

175 FORMAT! 4X, 1P6E16.6 ) 

176 FORMAT ( 1 X , I 3 , T 1 20 , ' A L PH A ' ) 

177 FORMAT (T120, 'BETA' ) 

210 FORMAT ( 'l',T5 ,'F(0,T0P)',T17,.'F( 1,T0P )',T29,'F(2,T0P)',T41, 

1 • F( 3, TOP ) ' ,T53, ' F(4,T0P ) ' ) 

12 FORMAT ( / 1P6E12.2) 

325 FORMAT ( / / , T,3 , ' ALT ' , T 1 2 , ' N ( HE ) 0 ' , T24 , ' N.( HE ) 1 ' , T36 , ' N ( HE ) 2 ' , 

1 T48, ' N ( H E ) 3 ' , T60, 'N(HE)4', T72, ' N ( HE ) 5 ' , T84 , ' HOR I Z . FLUX 

326 FORMAT (//,T3, 'ALT', T10, 

1 • HE ( 0 DEG) ' ,T22, 'HE(90 DEG ) ' , T 35 , ' HE ( 1 BOD EG ) ' , T 53,'RH(0)' 

1 T65 , ' RH ( 90 ) ' , T77, 'RH( 180) ' ,T95, 'FLUX(O) ' ,T107, ’FLUX(90) ' , 


OOOo 100U 

000U20O0 

0000 3 000 
00004000 
00005000 
00006000 
00007000 
00008000 
00009000 
00010000 
00011000 
00012000 
00013000 
00014000 
00015000 
00016000 
00017000 
00018000 
00019000 
00020000 
00021000 
00022000 
00023000 
00024000 
00025000 
00026000, 
00027000 
00028000 
00029000 
00030000 
00031000 
00032000 
00033000 
00034000 
00035000 
00036000 
00037000 
00038000 
00039000 
00040000 
00041000 
00042000 
00043000 
00044000 
00045000 
00046000 
00047000 
00048000 
00049000 
00050000 
00051000 
00052000 
00053000 
00054000 
00055000 
/ >00056000 
00057000 
00058000 
00059000 
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1 T1 19, 'FLUX! 180 ) 1 / ) 

1550 FORMAT! 1 X , -5P 1 F 5 . 0 , 1 P7E 1 2 . 2 ) 

151 FORMAT ( 1X.-5P1F5.0, 1P3E12.2, 6X , 1 p 36 12 . 2 , 6X , 1 P3E 1 2 . 2 ) 

400 FORMAT! //,T3, ' LAT 1 , T12» 'HE(120>', T24, 'HE<300)', 

1 T 36 * 1 HE ! 500 ) 1 / ) 

450 FORMAT (IX, 15, 1P3E12.2) 

WRITE! 6,5) 

C MUR IS THE NUMBER GF GRID POINTS ********************************* 
C NDIMl IS THE NUMBER OF HARMONICS OR DIMENSION OF THE MATRICES ***** 
C NT IS THE NUMBER OF TIME STEPS DESIRED *************************** 
C DELT IS THE VALUE UF THE TIME STEP ******************************* 
C NV IS ONE IN NORMAL TIME DEPENDENT C ALCUL AT I ON : I T IS SET NOT 
C EQUAL TO ONE IF STEADY STATE RESULT IS DESIRED 

C NWR IS ZERO IF NO PRINTOUT OF THE MATRICES IS DESIRED ************* 
C NF IS SET TO ZERO IF THE HORIZONTAL FLUX IN UPPER BOUNDARY IS 
C NOT DESIRED 

C NEDC IS SET TO ZERO IF EDC IS DESIRED TO BE CONSTANT 

DATA NT/ 1/, ND 1 M/6/ , NDR/420/, NV/O/, NWR/O/ 

DATA NF/1/, NEDC/ 1/ 

DELT=7.2D03*24.D0 
NP2=NDR+2 
NP 1 =NDR+ 1 
DO 225 M=1,NDIK 
PM ( M ) =0 . DO 
DPN(M)=0.D0 
DO 225 N = 1 , N D I M 
DO 22 5 L = 1 , ND I M 

CALL 0ATR(-1.00, l.DO, 1 .D-2,200, PLP, ALMN! L ,M,N ) , IER, AUX ) 

22 5 CALL OATR! -.99D0, .99D0, .0100,20, PL L ,BLMN( L,M,N ) , I ER, AUX ) 
BLMN(2,3,3)=0.D0 
BLMNl 2, 1 ,3 ) =0.000 
BLMN (2,2,2) = 0. 

BLMN! 2,3,1) = 0. 

PN( i)=l.D0 
PN< 3 ) =-0 . 5D0 
PN( 51=0.37500 
DPN ( 2 )=-l .DO 
DPN( 4 ) = 1 . 5D0 

DPN(6)=-1. 87500 
EDCM = 2.5D06 
ALT E 1 = 1.1D07 
ALT E2 = 1 . 3007 
DO 600 1=1, NP2 

R I = I 

ALT=ALTO+RI*DR 
HEIGHT! I ) =ALT 

EDC I ( I ) = EDCM *DEXP( 1 .10000D-06*! ALT- ALTE 1 ) )' 

IF (ALT.GE. ALTE 1 ) E DC I ( X ) = EDCM 

IF ( ALT.GE. ALTE2 ) EDC I ( I ) =EDCM*DEXP ( 1 . 10D-06* ( ALTE2-ALT ) ) 

IF (NEDC.EQ.O) EDC I ( I ) = 4.U06 
SCHT I ( I ) = DHH ( ALT ) 

SCHT A! I ) = SCHT! ALT, MBAR! ALT) ) 

DIF! I ) = D I FC ( ALT ) 

DENS! I ) = DENT(ALT) 

VLN = DENS! II 

CALL VLLIVL, BL, I, VLN, NDIM, 2, NV ) 

VLI(I) = VL ( 2 ) 

BLI(I) = BL ( 2 ) 

SHI(I) = SH( ALT ) 

SCIIII = SC ( ALT ) 


00060000 
00061000 
0006 2 QUO 
00063000 
00064000 
00065000 
00066000 
00067000 
00068000 
00069000 
00070000 
00071000 
00072000 
00073000 
00074000 
00075000 
00076000 
00077000 
00078000 
00079000 
00080000 
00081000 
00082000 
00083000 
00084000 
00085000 
000 8 60 00 
00087000 
00068000 
00089000 
00090000 
00091000 
00092000 
00093000 
00094000 
0UU95O00 
00096000 
00097000 
00098000 
00099000 
00100000 
00101000 
00102000 
UO 1 03000 
00104000 
00105000 
00106000 
00107000 
00108000 
00109000 
00110000 
00111000 
00112000 
00113000 
00114000 
00115000 
00116000 
00117000 
00118000 
00119000 


128 



ODH ( ALT ) 


DDH I ( I ) 

600 CUNTFNUE 

DO 501 13=1, NDIM 

DM ( I 3 ) =0. DO 
DO 501 14=1, NDK 


00120000 

00121000 

00122000 

00123000 

00124000 


BE ( 13, 1.4 ) = 0 • DO 
501 CONTINUE 

C BEGI N T I ME LOOP ^ — — — — — 

DO 503 I T= 1 , NT 

C SET BETA (BE) AT UPPER BOUNDARY TO ZERO 
DO 504 15 = 1 , NDIM 

B E ( I5,NDR)=0. DO 
504 CONTINUE 

C CALCULATE ALPHA ( AL ) AT UPPER BOUNDARY s*************************** 
WRITE (6,530) IT 

530 FORMAT (110X, 15H TIME STEP N0.=,I5) 

IP ( IT.GT.2) GO TO 550 
EDC=EDC I ( NDR ) 

CALL CALCltFM, BE ( 1 , NDR ) , AL ( 1 , 1 , NDR ) , NDR , ND I M , I T , NV , NF ) 

550 CONTINUE 
NDR 1=NDR-1 
DO 200 12=1, NDR] 

1= NP1 -1.2 


00125000 

00126000 

00127000 

00128000 

00129000 

00130000 

00131000 

00132000 

00133000 

00134000 

00135000 

00136000 

00137000 

00138000 

00139000 

00140000 

00141000 

00142000 


C 


R I = I 

ALT = ALTO + R I * DR 
DTSH = SCHTI ( I ) 

DTSC = SCHT A ( I ) 

DTDC = D I F ( I ) 

DTDD = ( DIF( 1 + 1 )-DIF( 1-1 ) )/( 2.D0*DR > 

EDC = EDC I ( I ) 

DTY = 1 . DO/ ( D 1 F ( I ) +EDC ) 

DTDH = DDH I ( I ) 

D1H= SCHT I ( I ) 

OSH = SHI ( I ) 

DSC = SCI ( I ) 

H = DENS ( I ) *2«*DR/(DENS( I — 1 ) -DENSII + 1)) 

R= RE + ALT 
V L ( 2 ) = VL I ( I ) 

BL ( 2 ) =BL I ( I ) 

I F ( IT. EQ. 1. AND.NV.EQ.l ) VL(2)=0.D0 

IF ( IT. GO. 1 .AND.NV.EQ.l ) BL(2)=0.D0 

GENERATE A, 8, C MATRICES ♦fc####**#*###*##!****’!'*#**#****###*##* 
DO 100 M = 1 .NDIM 
DO 100 N = l.NDIM 
SUM L = 0. 

SUMLB = 0. 

DO 11 L = l.NDIM 

10 SUML = SUML + VL ( L ) * ALMN ( L , N , M ) 

11 SUMLB = SUMLB + VL ( L ) * ALMN ( L , N , M ) / H +BL ( L ) * 

1 BLMN1 L , N »M ) / R 


00143000 
00144000 
00145000 
00146000 
00147000 
00148000 
00149000 
00150000 
00151000 
00152000 
00153000 
00154000 
00155000 
00156000 
00157000 
00 1 58000 
00159000 
00160000 
00161000 
00162000 
00163000 
00164000 
00165000 
00166000 
00167000 
00168000 
00169000 


RM = M-l 
C M N = 0 • 

I F < M.EO.N )CMN=-2.*RM*( Rh+1 . )/ ( 2.*RM + 1 . ) 
DE = 0. 

IF(M.EO.N) DE = 1 . 

RN=N— 1 

X = ( 2.4RN+1. )/2. 

BRA = (DTDD + DTDC* DSH + EDC / H ) V DE 


00170000 

00171000 

00172000 

00173000 

00174000 

00175000 

00176000 

00177000 


At N, M 
B ( N, M ) 


) = DE / (DR*DR) -(1. / (2.#DR)># ( BR A-X* SUML ) * DTY 00178000 

-2. * DE/I0R4DR) + DTY * ((DTDH + EDC * DTH ) * DE + DTDC 00179000 
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1 * CMN*X /( R* R) - SUMLB#X) 

C(N, M ) = DE / (DR*DR) + (1. / !2.*DR)) * uTY* ( BRA-SUML*X ) 
BP!N,M)= B ( N » l-i ) + DE* ( 2.D0/()ELT)*DTY 

I F ( IT.EQ. 1 ) GO TO 100 

B(N,M> = BIN, r-i ) -DE* ( 2 . DO/OEL T ) *1)T Y 

100 CONTINUE 

IF (NWR.EO.O), GO TO 520 
WRITE (6,251) 

WRITE! 6,250) A 
WRITE (6,252) 

WRITE (6,250) B 
WRITE! 6,253) 

WRITE! 6,250) C 
WRITE! 6, 176) I 

WR I TE ( 6 , 175 ) ( ( ALIN, M, I ) ,N=1,NDIM) ,K = 1,ND1M) 

520 CONTINUE 

C CALCULATE ALPHA ( AL ) FOR GRID PT. ( I — 1 ) 

C AFTER THE SECOND TIME STEP (IT.GT.2) ALPHA ( AL ) NEEIJ NOT BE RECAL- 


********** 


IF (IT.GT.2) 


GO TO 570 


CALL MINV! AL! 1, 1, I ) ,NDIM, D, LL , MV ) 

CALL GMPRD ( C , AL! 1,1,1), PRO, ND I M , ND I .M , ND 1 M ) 

CALL GMADD! PRO, B, XDEN, NO IN, ND Ii-'i ) 

IF (NWR.EO.O) GO TO 522 
WRITE! 6, 176) I 

WRITE! 6, 175 ) ( ( AL! N,ri, 1 ) ,N=1 ,NDIM) ,H = 1 ,NDIM) 

WRITE! 6, 51 ) 

WRITE (6,50) XDEN 
522 CONTINUE 

CALL SMPY! XDEN, -l.DO, DEN, ND I M , ND I M , 0) 

CALL SMPY! A, l.DO, UDEN, IMDIM, MOIM, 0) 

CALL MINV! A , l\'D IK , D, LL, MV) 

CALL GMPRD! A, DEN, A L ( 1 , 1,1-1), NO I M , NO IM , ND I M ) 

C CALCULATE D( DM) , ( EQUATION 4), «««*«*«?«» »** V V V V V v ¥ »,» V V ¥ v 3 S» V '1' 

570 CONTINUE 

C IF FIRST TIME STEP ( I T .EO . 1 ) ,SK I P CALCULATION OF 0** '■* V 3 i J ¥ t ¥ ¥ 3 £ ; r 'I 1 3 r 3 I s ¥ 3 1' 
I F ( IT.EO, 1 ) GO TO 510 

C AT THIS pt.,in THE FIRST and SECOND TIME STEPS (IT.LE.2),A Has 
C ALREADY BEEN INVERTED ******** »**»♦*«** *** **** * **** ** *****:* V s;c 

C AFTER THE SECOND TIME STEP, A HAS NOT BEEN INVERTED**************** 
C AT THIS POINT, UDEN IS THE SAME AS THE UNI VERTED A MATRIX ******** 
IF (IT.GT.2) GO TO 610 

CALL GMPRD ( UDEN , F ( 1 , I - 1 ) , DM, NDIM, NDIM, 1) 

GO TO 620 
610 CONTINUE 

CALL GMPRD (A, F ( 1 , I — 1 ) , DM, NDIM, NDIM, 1) 

620 CONTINUE 

CALL GMPRD ( BP , F ( 1 , 1 ) , DM1, NDIM, NDIM, 1) 

CALL GMPRD! C, F( 1,1+1), DM2, NDIM, NDIM, 1) 

CALL GMADD (DM, DM1, DK3 , NDIM, 1) 

CALL GMADD (DM3, DM2, DM1, NDIM, 1) 


F ( 1,1-1) , DM 


F ( 1 , 1 ) , DM1, NDIM, 

F( 1,1+1), DM2, NDIM 
DM1, DM3, NDIM, 1) 
DM2, DM1, NDIM, 1) 


CALL SMPY (DM1, -l.DO, DM, NDIM,1, O) 

C CALULATE BETA (HE) AT GRID PT. ( I — 1 > ****************************** 
C ALPHA! I — 1 ) HAS NOT BEEN INVERTED YET, ALPHA! I) HAS BEEN INVERT El)** 


CALL GMPRD (C, AL! 1,1,1), ALP, NDIM, NDIM, NOI 

CALL SMPY (ALP, -l.DO, UDEN, NDIM, NDIM, 0) 

CALL GMPRD (UDEN, BE! 1,1), NUN, NDIM, NDIM, 1) 

CALL GMSUB (NUM, DM, ALF, NDIM, 1) 

CALL SMPY (ALF, -l.DO, PROB, NDIM, 1, 0) 

IN THE FIRST AND SECOND TIME ST E P S ( I T . L E . 2 ) A HAS ALREADY BEEN 
INVERTED: SO SKIP TO 630 ********************************* 
IF ( IT.LE.2 ) GO TO 630 

CALL MINV (A, NDIM, D, LL, MV) 

630 CONTINUE 

CALL GMPRD (A, PK06, 6E(1,I-1), NDIM, NDIM, 1) 

510 CONTINUE 
200 CONTINUE 


NDIM, NDIM, NDIM) 
NDIM, 0) 

NO I M , 1 ) 


'I s ¥ 3 ! C 'i' 5 i ; 'I* 3 0 3 i S 3 i' 3 I' 3 i 5 ¥ 3 i : 


00180000 
00 181000 
00 182000 
001 83000 
00 184-000 
ooi a 5 ooo 
00186000 
00187000 
001.88000 
00189000 
00190000 
00191000 
00192000 
00193000 
00 1 94-000 
0019 5 0 0 0 
00196000 
00197000 
00198000 
0 0 1 9 9 0 U 0 
00200000 
00201000 
00202000 
00203000 
00204-000 
00205000 
00206000 
00207000 
00208000 
00209000 
0021.0000 
00211000 
00212000 
00213000 
00214000 
00215000 
00216000 
0021 ( 000 
00 2. 1 8000 
00219000 
00220000 
00221000 
00222000 


DM, 

NDIM, 

ND I M , 

1 ) 

00223000 

002241)00 

00225000 

DM, 

NDIM, 

N DIM, 

1 ) 

00226000 

00227000 

DM1, 

NDIM, 

N D 1 M , 

1 ) 

00228000 
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C THIS IS THE END OF THE INTERMEDIATE STEPS : THE DENSITIES WILL NEXT BE 
C CALCULATED 

C SET BOUNDARY CONDITION AT LOWER BOUNDARY ***¥**S4SS*******Y#*Y#*4 
F ( 1,1) = 1.663D09 

DO 497 11 = 2* NO I M 

F( 1 1, 1 >=0.00 


00248000 

00249000 

00250000 

00251000 

00252000 

00253000 


497 CONTINUE 

C INVERT ALPHA (AL) AT LOWER BOUNDARY; THE REST OF THE AL'S HAVE 
C ALREADY BEEN INVERTED ************* Jit****# #44*’? *>? sc#*#*:#**#!**# 

C AFTER THE SECOND TIME STEP, THIS NEED NOT BE DONE 
IF ( IT.GT.2) GO TO 580 

CALL M I N V ( A L ( 1 * 1 * 1 ) * MO I M * 0, LL, MV) 

580 CONTINUE 

DO 300 J = 1 * N D k 
K = J 

CALL GMSUB ( F ( 1 , K ) , BE(1,K), DM4, NDIM, 1) 

CALL GNPRD < A L ( 1,1 ,K) , DM4 , F ( 1 , K + l ) , NDI M . ND I M, 1 ) 

DH 1 ( K ) = DHEIO.ODO, K) 

DH2 ( K ) = DHE( 1.5708 DO, K) 

DH3 ( K ) = DHE ( 3 . 1416D0, K) 

ZETKK) = DH1 ( K ) / ( DENS ( K ) + DH 1 ( K ) ) 

Z E T 2 ( K ) = DH2 ( K ) / ( DENS ( K ) + DH2 ( K ) ) 

ZET3IK) = DH3 ( K ) / ( DENS ( K ) + DH3 ( K ) ) 

300 CONTINUE 

HEIGHT! 1 )=81 .E05 
WRITE! 6,325 ) 

B I GPH I =0 .DO 
DO 125 1= 1 , N D R 
R = RE + HE IGHT! I ) 

PHI=0.D0 

PH I 1 = 0.00 

DO 128 1 9=1, NDIM 

PH I = PHI + F( 19, I )*OPN( 19) 

PHI1= PH 1 1 + F ( 19, I )*PN( 19) 

128 CONTINUE 

PHI= — PHI *0 I F ( I ) / K 
PH I 1 = -PHI l*t)LI ( I )*0.5D0 

PH I = PH I + PH I I 

BIGPHI=BIGPHI + PHI*1 .1)05 

125 WRITE! 6 , 150) HE I GHT ( I ) , ( F ( J , I 1 , J = 1,NDIM) , PHI 
WRITE! 6,326) 

DO 126 1=1, NDR 

IF (I .GE.NDR1 ) GO TO 127 

IF ( I.EO. 1 ) GO TO 127 

FLU X 1 = - ( DIF! I ) +EDCI! 1 ) ) # ( 0H1 ( I + 11-DH1! I — 1 ) )/(2.00*DR) 
1-DHKII * ( S H I ( I ) *D I F ( I ) +SCI(I) *EDCI(I) )+DHl(I) *VLI(I) 

FLUX2= -( DIF(I) +EDCI! I ) )*(DH2( I + D-DH2! 1-1) )/(2.D0*DR) 

1 -DH2II) * ( S H I ( I ) -DIF! I ) +SCIIII *EDCI(D) 

FLUX3= -( DIF! I ) +EDCI! I) ) * ( DH3 ( I + 1 ) -DH3 ( 1-1) )/(2.D0*DR) 

1-DH3 ( I ) * ( S H I ( I ) *D I F ( I ) +SCIIII *EDCI(I) )-DH3(I) *VLI(I) 

GO TO 126 
127 FLUX1=0.00 
FLU X2 = 0 . DO 
FLUX3=0.D0 

126 WRITE (6, 151) HE I GHT ( I ) , 

1 DH1 ( I ) , DH2! I ) ,DH3< I ) , ZETKI), ZET2(I), ZE T 3 ( I ) ,FLUX1, 


00254000 

00255000 

00256000 

00257000 

00258000 

00259000 

00260000 

00261000 

00262000 

00263000 

00264000 

00265000 

00266000 

00267000 

00268000 

00269000 

00270000 

00271000 

00272000 

00273000 

00274000 

00275000 

00276000 

00277000 

00278000 

00279000 

00280000 

00261000 

00282000 

00283000 

00284000 

00285000 

00286000 

00267000 

00288000 

00289000 

00290000 

00291000 

00292000 

00293000 

00294000 

00295000 

00296000 

00297000 

00298000 

00299000 

00300U00 

00301000 

00302000 

00303000 
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1 FLUX2.FLUX3 

WRITE! 6,400) 

DO 350 11=1,19 
R I = I I 

DHL1 = DHE( ( RI-1 .L>0)*0. 17453300,40) 

IF ( II.EO.l) R 1 = DHL 1 
IF ( 1 I .EO. 19) R2=DHL 1 

DHL2 = D HE { (RI-1. D0)*0. 17453300,220) 

IF (II .EO. 1 ) R3 = DHL2 

IF (II.E0.19) R4=DHL2 

DHL3 = DHE( ( RI-1 .DO )40.174533D0,420) 

IF (II.EO.l) R5=DHL3 
IF (II.E0.19) R 6 = 1) FI L 3 
LA = 10*< I I -1 ) 

350 WRITEI6, 450) LA ,DHL1 , DHL2,OHL3 
R1=R2/R1 
R2=R4/R3 
R 3 = R 6/R5 
WRITE (6,130) 

130 FORMAT! //T12, 'RATIOS OF POLE OENS. ' ,T60, ' INTEG. FLUX AT EO.'/) 
WRITE (6,131) Rl, R2,R3, BIGPHI 

131 FORMAT ( 6X , 1 P 3D 1 2 . 2 , T60 , 1 P 1 D1 4. 4// ) 

503 CONTINUE 

WRITE (6,650) 

WRITE (6,652) 

651 FORMAT (IX, -5P1F5.0, 1 P6D 1 4 . 2 , - 5P2F 1 4 . 4 , 1 P 1 U 1 2 . 2 ) 

650 FORMAT! ////,T3, 'ALT(KM) 1 ,T15, 1 EDC * , T29 , ' D I FC ' , T43 , 

1 ' VL/ ' ,T57, • HL/' , T71 , ’ VL ' , T85, 1 BL ' ., T99 , ' 1 /SC ' , T113, 

1 1 1 /SH ' , T 1 2 5 , 1 X/N ' ) 

652 FORMAT (T13, • ( CM*CM/SEC ) ' , T27 , ' ( CM*CM/SEC ) ' , T41 , 

1 1 ( SH4DIFC )' ,T55, '( SH*DIFC )' ,T69, • (CM/SEC )' ,T83, ' (CM/SEC) ' , 

1 T99 , '(KM)', T 1 1 3 , '(KM)'/) 

X8=DENS< 2201*100. DO 
DO 640 17=1, NDR 

ALT = HE I GHT ( 17) 

EDC=EDC I (17) 

OTOC=D I F ( 17) 

VLHD=VL I ( I7)/(SHI( I7)4DTDC> 

BLHD = BL I ( I 7 ) / ( SH I ( I7)*DTDC) 

R I NV 1 = 1 . DO/ SC I ( 17) 

R1MV2=1. DO/SHI ( 17) 

X9=X8/DENS( 17) 

WRITE (6,651) ALT, EDC, DTDC , VLHD , 8LHD , V L I ( 1 7 ) , hi L I ( I 7 ) 
1 , R I NV 1 , RINV2 ,X9 

640 CONTINUE 
STOP 
END 

BLOCK DATA 

IMPLICIT R E A L* 8 ( A-H,0-Z) 

COMMON /MOD/ TEMP ( 41 ) , MM ( 41 ) , DENN2 ( 4 1 ) ,DEN02(41 ) ,DEN01 ( 41 ) , 

A DENHEI 41 ) ,DENA( 41 ) ,MN2 ,M02,M01,MHE ,MA 
REAL*8 

1 MASS, MM, M N 2 , M02 , M01, MHE , MBAR 

DATA TEMP / 3* 1 86 . 0 , 5* 1 8 5 . 9 , 2* 1 8 5. 8 , 1 B8 . 4 , 1 90 . 9 , 1 9 3 . 5 , 1 9 5 . 9 , 

A 198.2,200.4,202.4,204.4,206.3,208.1,212.2,215.7,220.0,224.6, 

B 229.0,233.4,237.9,242.3,246.8,251.1,261.6,271.9,282.3,292.7, 

C 302. 9, 313. 1,323. 6, 334. 0,344. 4, 355.0/, MM /4*28 . 96 , 5*28 . 9 5 , 

D 2428.94,28.92,28.89,28.87,28.83,28.78,28.70,28.61,28.52,28.42, 


00304000 

00305000 

00306000 

00307000 

00308000 

00309000 

00310000 

00311000 

00312000 

00313000 

00314000 

00315000 

00316000 

00317000 

00318000 

00319000 

00320000 

00321000 

00322000 

00323000 

00324000 

00325000 

00326000 

00327000 

00328000 

00329000 

00330000 

00331000 

00332000 

00333000 

00334000 

00335000 

00336000 

00337000 

00338000 

00339000 

00340000 

00341000 

00342000 

00343000 

00344000 

00345000 

00346000 

00347000 

00348000 

00349000 

00350000 

00351000 

00352000 

00353000 

00354000 

00355000 

00356000 

00357000 

OU 3 5 8000 

00359000 

00360000 

00361000 
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E 28.30,28.18,28.02,27.99,27.92 
F 27.66,27.57,27.49,27.41,27.34 
G/,DENN2 / 


27.87,27.82,27.78,27.74,27.71, 
27.26,27.19,27.13,27.08,27.04,27. 
2.478E 14, 2 .072E 14.1.733E 14,1.449E 14, 


H 

1.212E 

14, 1.014E 

14,8.480E 

13, 7.095E 

13, 5.934E 

13,4.965E 

13, 

I 

4 . 1 0 3 E 

13,3. 544E 

13, 2. 83 IE 

13,2 .349E 

13, 1.947E 

13, 1.626E 

13, 

J 

1.362E 

13, 1. 146E 

13, 9.673E 

12,8. 178E 

12,6.8176 

12, 5.704E 

12, 

K 

4.804E 

12 ,4.060E 

12,3 .453E 

12,2 .950E 

12,2.529E 

1 2 , 2. 174E 

12, 

L 

1.875E 

12,1. 620E 

12, 1 .365E 

12,1 . 164E 

12,9. 983E 

11,8. 606E 

11, 

Mi 

7.460E 

11 ,6.513E 

11 ,5.723E 

11,5.0576 

11,4.4786 

1 1 ,4.008E 

11 / 

DATA DEN02 / 

6. 

649E 13,5. 

559E 13,4. 

648E 13,3. 

8886 

A 

3.251E 

13,2. 721E 

13,2.275E 

13, 1.906E 

1 3 , 1 . 598 E 

13, 1.332E 

13, 

B 

1. 101E 

13,9. 188E 

12,7.3616 

12,6. 146E 

12,5.1906 

1 2 , 4. 296E 

12, 

C 

3.553E 

12 ,2 .936E 

12 ,2.423E 

12,1.9946 

12 , 1.644c 

12, 1.359E 

12, 

D 

1. 131E 

12,9. 443E 

11,7.9326 

1 1,6.6936 

11,5.6656 

11,4. 809E 

11, 

E 

4.093E 

11,3 .492E 

11,2.9036 

1 1 ,.2 .443E 

11 ,2.066E 

11, 1.757E 

11, 

F 

1.501E 

11, 1.292E 

11,1.1196 

11,9.7 44E 

10,8. 501E 

10, 7.49 5 E 

10/ 

G 

, D E N 0 1 

/ 

8.700E 10,8.9306 10.9.210E 10,9.5006 10, 

H 

9.800E 

10, 1.015E 

■ 1 1 , 1 • 0-55E 

11,1.1056 

11,1.1656 

11, 1.250E 

11, 

I 

1.420E 

11,1 .680E 

11 ,2.060E 

11,2.6606 

11,3.4106 

11,4. 100E 

10, 

J 

4.515E 

11,4. 800E 

11,4.935E 

11, 5. OOOE 

11,4.9456 

1 1 , 4 . 7 6.0 E 

11, 

K 

4.425E 

11 ,4.050E 

11,3.6106 

11, 3.21 OE 

11,2.8356 

11,2.5106 

11, 

L 

2.230E 

11,2. OOOE 

11, 1.812E 

11, 1.642E 

1 1 , 1.487E 

11,1.3476 

11, 

Ml 

1.235E 

11,1.125E 

11,1.0206 

11,9. 250E 

10,8.400E 

10,7. 600 E 

10 / 

DATA DENHE / 

1 . 

663E 09,1. 

39 IE 09,1. 

163E 09,9. 

725E 

A 

8. 132E 

08,6. 807 E 

08,5. 69 1 E 

08,4.7616 

08,3.9826 

08,3.3326 

08, 

B 

2.753E 

08,2. 282E 

08, 1 .886E 

06,1. 568E 

08,1. 306E 

08, 1 .091E 

08, 

C 

9. 144E 

07 , 7 .692E 

07,6 .49 2 E 

07, 5.492E 

07 , 4 . 57 5E 

07 , 4.42 1 E 

07, 

0 

4.275E 

07,4. 1 3 8 E 

07,4.0086. 

07,3. 886E 

07,3.7.68E 

07,3.6566 

07, 

E 

3.551E 

07.3.450E 

07,3.3046 

07,3.1716 

07,3.048E 

07,2.934E 

07, 

F 

2.829E 

07, 2. 73 IE 

07,2. 639E 

07, 2. 554 E 

07,2.4746 

07 , 2 . 400E 

07/ 

G 

, D EN A 

/ 

2.965E 

12.2.479E 

12,2.0736 

12, 1.734E 

12, 

H 

1.450E 

12, 1.213E 

12, 1.014E 

12, 8.486E 

11,7.0986 

1 1 , 5 . 939E 

11, 

I 

4.868E 

11 ,4.046E 

11,3.363E 

11,2.7956 

11 ,2.329E 

11, 1.945E 

11, 

J 

1.630E 

11, 1.371E 

11,1.1576 

11,9.8006 

10,8. 1 54E 

10,6.8236 

10, 

K 

5.746E 

10, 4. 85 7 E 

10,4. 130E 

10,3.5286 

10,3.0256 

10,2.6016 

10, 

L 

2.242E 

10, 1,9386 

10, 1 .633E 

10, 1 . 393E 

10, 1.1946 

10, 1 .029E 

10, 

Ml 

8.923E 

09 , 7 . 791 E 

09,6.8466 

09,6.0496 

09,5.3576 

09,4.7956 

09 / 


END 

BLOCK DATA 
IMPLICIT 


REAL*8(A-H,0-Z) 


COMMDN/CALC/DR, CS, RE, ALTO,' TINE, EDC , TLB , ZLB, M.N2, NO 2, 

1 M 0 1 , MHE’ 

REALMS 

1 MN2, M02, M01, MHE 

DATA DR/ 1 . 00 E 05/ ,CS/ 2 . 12E- 1 5/ , RE /6356 . 77E 05 / , ALTO/80 . OOE 05/, 

1 TIIMF/1100./,EDC/4.00E 06/ , TLB/ 355. / , ZL B/ 120 . OOE 05/ , MN2/ 28. / , 

2 M 02/32./, M01/16./ ,MHE/4./ 

END 


FUNCTION MBAR ( ALT ) 

IMPLICIT RE A L*8 ( A-H,0-Z ) 

R E A L * 8 

1 NASS, MM, MN2 , K02,'M01, MHE, MBAR 

MEAN MASS FROM 80 KM TO TOP 

COMMON/CALC/DR, CS, RE, ALTO, TINF, EUC, TLB, ZLB, MN2, M02, 

1 ,M 01, MHE 

COMM ON /MOD/TEMP ( 41 ) , MMI41), DENN2I41), DEN02I41), DEN01 ( 41 ) , 
1 DENHE ( 41 ) , DEIMA141) 

COMMON/INDEX/L , M, N, I, LI, Ml, N1 
REAL-8 

1 MASS, MN2 , M02 , M01, MHE, MBAR 


00362000 
0100363000 
00364000 
00365000 
00366000 
00367000 
00368000 
00369000 
00370000 
, 00371000 
00372000 
00373000 
00374000 
00375000 
00376000 
00377000 
00378000 
00379000 
00380000 
00381000 
00382000 
00383000 
00384000 
, 00385000 
00386000 
00387000 
00388000 
00389000 
00390000 
00391000 
00392000 
00393000 
00394000 
00395000 
00396000 
00397000 
00398000 
00399000 
00400000 
00401000 
004U2000 
00403000 
00404000 
00405000 
00406000 
00407000 
00408000 
00409000 
00410000 
00411000 
00412000 
00413000 
00414000 
00415000 
00416000 
00417000 
00418.000 
00419000 
00420000 
00421000 


133 



IF (ALT - ZLB ) 25, 25, 45 
25 MBAR = M M ( I ) 

GO TO 50 

45 MBAR = ( DN2 ( ALT ) * MM2 + D02 ( ALT ) * M02 + DOl(ALT) * hUl + 

X DHE(ALT) 4 MHE ) / ( DN2 ( ALT ) + D02 ( ALT ) + DOKALT) + DHE ( ALT ) ) 
50 RETURN 
END 

SUBROUTINE VLL(VL, BL, I I , VLN ,ND I M , I T , N V ) 

IMPLICIT REAL*8(A-H,0-Z) 

DIMENSION V L < 3) , BL ( 3 ) 

COMMON/ INDEX/ L ,M,N , I, LI, Ml, NX 

COMMON/CALC/DR, CS, RE, ALTO, TINF, EDC , TLB, ZLB, MN2, IT02, 

1 M01 , MHE 

REALMS 

1 MASS, MM, MN2 , M02 , HOI , MHE, MBAR 

DT ( ALT ) = (TIALT + DR) - T ( ALT -DR ) ) / ( 2 .-DR ) 

TDT( ALT ) = DT( ALT ) / T ( ALT ) 

SCHT ( ALT , MASS ) = 8. 31E07*T ( ALT )*( (RE+ALT )**2 ) /( MASS-980. 665* ( R E 
1 **2 ) ) 

SC ( ALT ) = TDT ( ALT ) + 1 ./ SCHT ( ALT , MBAR < ALT ) ) 

DENT ( ALT ) = DN2 ( ALT ) + D02 (ALT) + DOl(ALT) 

DHH(ALT) = (SC (ALT + DR) SC 

1 ( ALT - DR ) ) / ( 2 • *DR ) 

R I = II 

ALT = ALTO + RI * OR 
D SC = DHH ( ALT ) 

SC1=SC( ALT ) 

R = RE + ALT 

H = DENT ( ALT )*2.*DR/( DENT ( ALT-DR) -DENT ( ALT + DR ) ) 

ZB = 80.D05 

ZT = 602. D05 
VW = 100. DO 

BETA = 1.8D-07 
DO 56 11=1, NDIM 
VL ( II) =0. 

BL ( I 1 ) =0 . DO 
56 CONTINUE 

IF ( IT.EO.l.AND.NV.EO.l ) GO TO 50 
IF (ALT - ZB) 10, 20, 20 
10 VL ( 2 ) = 0. 

BL ( 2 ) =0 . DO 
GO TO 50 

20 IF (ALT - ZT) 25, 35, 35' 

25 DX=< ALT-200. D05 ) 

ALN= BETA *DX 

VL ( 2 ) = ( VW/2 . ) * ( 1 . +DERF ( ALN ) ) 

X2=(VW/1. 77245) * BETA *DEXP ( -ALN*ALN). 

BL ( 2 ) =R* ( X2 -VL( 2 )*SC1 + 2.*VL ( 2 )/R ) 

GO TO 50 

35 VL<2) = VW 

BL ( 2 ) =V L ( 2 ) *( 2 .- R*SC 1 ) 

50 RETURN 
END 

SUBROUTINE CALC 1 (FM, BET, AL PH A , NDR , ND I M , I T , NV , NF ) 

IMPLICIT REAL*8( A-H ,0-Z ) 

REAL*8 

1 MASS, MM, MN2 , M02 , M01, MHE, MBAR, NUM 


00422000 

00423000 

00424000 

00425000 

00426000 

00427000 

00428000 

00424000 

00430000 

00431000 

00432000 

00433000 

00434000 

00435000 

00436000 

00437000 

00438000 

00439000 

00440000 

00441000 

00442000 

0044300(i 

00444000 

00445000 

00446000 

00447000 

00448000 

00449000 

00450fi00 

00451000 

00452000 

00453000 

00454000 

00455000 

00456000 

00457000 

00458000 

00459000 

00460000 

00461000 

00462000 

00463000 

00464000 

00465000 

00466000 

00467000 

00468000 

00469000 

00470000 

00471000 

00472000 

00473000 

00474000 

00475000 

0047600.0 

00477000 

00478000 
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CflMMON/CALC/DR , CS, KE, ALTO, TIME, EOC, TLB, ZLB, MN2, M02, 

1 KOI, H E 

COMMON/ INDEX/L, M, N , I, LI, Ml, N1 
COMMON / CO E F / 6LMN(6,6,6), BLMN( 6,6,6) 

DIMENSION V L ( 6 ) , HL(6) 

DIMENSION U AMM (6,6), UDENF(6,6) 

DIMENSION X AMM ( 6,6) ,XDENF( 6,6 ) 

DIMENSION L L A ( 6 ) , KM A ( 6 ) , LLN(A), MMN ( 6 ) 

DIMENSION A M K ( 6,6 ) , BMi"i( 6, 6 ) , ABM ( 6 , 6 ) , ADM { 6 ) , NUriF ( 6 ) , DM ( 6 ) , BE T ( 6 ) , 
1 AL(6,6,1),DENF(6,6),FM(6), ALPHA(6,6) 

REAL-8 NUMF 

DTI ALT) = (TIALT + DR) - TIALT - OR ) ) / ( 2 ,-IJK ) 

TDT(ALT) = DTI ALT ) / T ( ALT ) 

SCHT I ALT, MASS )= 8.31E07*T ( ALT )*( I RE + ALT )--2 ) / I MASS-9 80 . 66 5- I KE 
1 ** 2 ) ) 

SH I ALT ) = (1. - .A ) -TDT I ALT ) + 1 . /SCHT I ALT , MHE ) 

SCI ALT >=TDT< ALT ) + 1 . DO/ SCHT I ALT , MBAR I ALT ) ) 

DENT! ALT )=DN2( ALT 1+0011 ALT 1+D02I ALT ) 

DIFCI ALT ) = ( 1.69D 19/DENT ( ALT >)*(( T < ALT1/273. 16)-- 0.691) 

FNDR=NDK 

ALT 1 =ALT O+FNDR-DR 
A L T 2 = A L T 0+ I F ND R+ 1 . DO ) *DR 
NP 1 =NDR+ 1 
D A = D I F C I ALT1 ) 

D B = D I F C ( ALT2) 

DC=1 .DO/ I DA+EDC ) 

DD= 1 . DO/ I DB+EDC ) 

SH 1 = SH I A LT 1 ) 

SH2=SH( ALT2) 

SC 1 = SC I ALT 1 ) 

SC2 = SC I ALT2 ) 

WRITE (6,701) DA, 08 , DC, DO 

WRITE (8,701) SHI, SH2 , SCI, SC 2 

701 FORMAT (D20.10, 020.10,020.10, 020.10) 

R 1 =Rt + ALT 1 

R2=RE+ALT2 

AM = — l.DO/OR +DA*SH1-DC/2.D0+EDC*SC1*DC/2.D0 
BM = l.DO/DK +DH«SH2 *00 / 2 . DO + EOC- SC2-OD/2 . 00 

WRITE (6,702) AM, fil-l 

WRITE (8,702) AM, BK 

702 FORMAT (5X, 4H AM= , 020. 10, 4H BM= , 020 . 10 ) 

DATA DM/ 3*0 . / 

V L N 1 = DENT ( ALT 1 ) 

V L N 2 = DENT ( ALT2 ) 

CALL V L L < VL , BL , NOR, VLN1 ,NDIM,IT,NV ) 

VTL=VL( 2 ) 

CALL VLl< VL,BL,NP1 , VLN2,NDIM, IT, MV ) 

VTT = VL ( 2 ) 

C WRITE (6,703) A LMN 

C WRITE (6,704) BLBN 

703 FORMAT (10X,6H ALMN= , 3020 . 10 ) 

704 FORMAT ( 5 X , 6H BL MM= , 3020 . 10 ) 

WRITE (6,707) VL 

WRITE (8, 707) BL 
707 FORMAT ( 5X,6H V L BL = , 3D20 . 10 ) 

WRITE (6,707) BL 
WRITE (8,707) VL 
C MAKE Ah A MATRIX = AMM ( N , M ) 

C Ml A K E 8 M A M A j K I X = 5 M M ( N , Ft ) 

IF (NF.EQ.O.AND. IT. £0.1 ) GO TU 710 


00479000 
00480000 
00481000 
00482000 
00483000 
00484000 
00485000 
00486000 
00487000 
00488000 
00489000 
00490000 
00491000 
00492000 
00493000 
00494000 
00495000 
00496000 
00497000 
00498000 
00499000 
00500000 
00501000 
00502000 
00503000 
00504000 
00505000 
00506000 
00507000 
00508000 
00509000 
00510000 
00511000 
00512000 
O0513000 
00514000 
00515000 
00516000 
00517000 
00518000 
00519000 
00520000 
00521000 
00522000 
00523000 
00524000 
00525000 
00526001) 
0O5270OU 
00528000 
00529000 
00530000 
0053 1000 
(MJ 532000 
00533000 
00534000 
00535000 
00536000 
00537000 
00538000 
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EPS1 = K1 / SCHT( ALTl.MHE ) 00539000 

FPS2 = K2 /SCHT ( A L T 2 » H H E ) 00540000 

Vfci Ak 1 = 1.D04* DSOkT (0.6192D04T ( 4 L T 1 ) ) 00541000 

VBAR2 = 1 .004 *DSORT ( 0.6192D0*T( ALT2 > ) 00542000 

HF1 = ( 1 . DO + ( H.4D0/EPS1 ) )4=VBAR1 /( EPS 1**2 > 00543000 

HF2 = (1.00 + ( 8.4D0/EPS2 ) )*VBAR2 /<EPS2**2) 00544000 


WRITE (6,902) HF1 
WRITE (8,902) HF1 
HF1=HF1*DC 
HF2=HF2-0D 
DO 11 M= l.NDIM 
DO 11 N = l.MDIK 
R i v i = M - 1 
X=RM*( Rl'i+l ) 

DELTA = 0. 

I F( N.EO.M DELTA = 2.D0/( 2 . DO*RN+ 1 . DO ) 

AMM ( N * M ) = ( AM + X4=HF 1 / 2 . DO ) *DELTA - ( DC/ 2 .DO ) 4=VTL*ALMN ( 2 , N , H ) 
11 BNh ( N »M ) =( BM+X4=HF2/2.DO)*DELTA H DD / 2 . DO ) *=VT T* ALMN ( 2 , N , M ) 
GO TO 720 
710 CONTINUE 

DO 10 M= 1 , N D I K 
DO 10 N = 1 , ND I M 
RM=K-1 
DELTA = 0. 

IF(M.EO.M) DELTA=2. 00/(2. DO*Rrt+l. DO) 

A M M ( N , h ) = AM 4= DELTA - ( DC/2 . DO > 4=VT L 4= ALMN ( 2 , N , M ) 

10 BMM ( N* M ) = BM * DELTA - ( DD/ 2 . DO ) 4=VTT 4=ALHN ( 2 , N , M ) 

720 CONTINUE 

WRITE (6,708) AMM 
WRITE (6,708) Bi-lK 


00545000 

00546000 

00547000 

00548000 

00549000 

00550000 

00551000 

00552000 

00553000 

00554000 

00555000 

00556000 

00557000 

00558000 

00559000 

00560000 

00561000 

00562000 

00563000 

00564000 

00565000 

00566000 

00567000 

00568000 


708 FORMAT (10X.8H AMMBMK= , 3020 . 10 ) 

WRITE (8,708) AMM 
WRITE (8,708) 6 Hi-1 

CALL MINVIAMM, NDIM, DETF , LLN, MMN ) 

902 FORMAT ( 5 X , 3H K=,1R1D13.5) 

CALL Gi'iPRD ( AMM , BMh , UDENF , NDIM, NDIM, IMDIM) 

CALL SMPY(UDENF,-1. DO, ALPHA, NDIM, NDIM, 0) 

WRITE (6,901) ALPHA 
901 FORMAT (5X, 7H ALPHA*, 3020. 10) 

WRITE (8,901) ALPHA 

RETURN 

END 

FUNCTION DN2( ALT ) 

IMPLICIT REAI 4=r( a-h,U-Z ) 

REAL w 8 

1 MASS, MM, MN2, M 02, Mill, MHE, MBAR 

M2 DENSITY FROM 80 KM TU TOP 

COM, MON/CALC/OR , CS, RE, ALTO, TINF, EDC , TLB, ZLB, MN2, M02, 

1 Mdl, MHE 

COMMON/f-iOD/TEMR ( 41 ) , MM(41), DENN2(41)» 0EN02141), DEN0K41), 
1 D E N H E ( 4 1 ) , DEN A ( 4 1 ) 

I = (ALT -ALTO) / DR + .5 
IF (ALT - ZLB) 25, 25, 45 
25 DIM 2 = DENN2 ( I ! 

GO TO 50 

45 DN2 = D J N 2 ( ALT ) 

50 RETURN 
END 

FUNCTION D02( ALT ) 

■IMPLICIT REAL4=8( A— H ,0-Z ) 


00569000 

00570000 

00571000 

00572000 

00573000 

00574000 

00575000 

00576000 

00577000 

00578000 

00579000 

<10580000 

ooooTooo 

00002000 

00003000 

00004000 

00005000 

00006000 

00007000 

00008000 

00009000 

00010000 

00011000 

00012000 

00013000 

00014000 

00015000 

00016000 

00017000 

00018000 
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R E A L * 8 





000 19000 

1 MASS, Mi-1, MN2, M02, M01, MHE, MHAk 





00020000 

02 DENSITY FROM 80 KM TO TOP 





00021000 

COMMON/CALC/DR, CS, RE, ALTO, TINF, EOC, 

TLH, 

ZLB, 

MN2, 

MO 2, 

00022000 

1 M01 , MHE 





00023000 

COMMON /MOD / T EM P ( 4 1 ) , MM (41), DENN2I41), 

DENU2 ( 41 ) , 

DENU1 (41 ) , 

00024000 

1 OENHEI 41 ) , DENAI41 ) 





00025000 

I = (ALT -ALTO) / UR + .5 





00026000 

IF (ALT - ZLB) 25, 25, 45 





00027000 

2 5 D02 = DEMO 2 ( I ) 





00028000 

GO TO 50 





00020000 

45 D02 = DJ02IALT) 





00030000 

50 RETURN 





00031000 

END 





00032000 

FUNCTION D01( ALT) 





00033000 

IMPLICIT REAL*8( A-H,0-Z) 





U0034000 

REALMS 





00035000 

1 MASS, MM, MM2, M02 , M01 , MHE , MHAk 





00036000 

01 DENSITY FROM 80 KM TO TOP 





00037000 

COMMON/CALC/DR, CS, RE, ALTO, TINF, EDC , 

TLH, 

ZLB, 

MN2, 

MO 2, 

00038000 

1 Mill, MHE 





00039000 

COMMON/HOD/TEMPI 41 ) , MM(41), DENN2I41), 

DENU2 ( 41 ) , 

DEN01 (41 ) , 

00040000 

1 DENHE ( 4 1 ) , DENAI41) 





00041000 

I = (ALT -ALTO) / DR + .5 





00042000 

IF ( ALT - ZLH ) 25, 25, 45 





00043000 

25 DO 1 = DEN01I I ) 





00044000 

GO TO 50 





00045000 

45 DO 1 = D J01 ( ALT ) 





00046000 

50 RETURN 





00047000 

END 





00048000 

FUNCTION DHE ( ALT ) 





00049000 

IMPLICIT REAL*8( A-h,0-Z ) 





000 500 UO 

REAL*8 





00051000 

1 MASS, MM, WN2, M 02, MG1 , MHE, MHAK 





00052000 

HE DENSITY FROM 80 KM TO TOP 





0U053000 

COMMON/CALC/DR, CS, RE, ALTO, TINF, EDC, 

TLB, 

ZLB, 

MN2, 

MO 2, 

00054000 

1 M 01, MHE 





00055000 

COMMON / MOD / T EMP ( 41 ) , MM(41), DENN2(41), 

DEN02 ( 41 ) , 

DENG 1(41 ) , 

00056000 

1 DENHE ( 4 1 ) , DENA ( 4 1 ) 





00057000 

I = (ALT -ALTO) / DR + .5 





00058000 

IF (ALT - ZLB) 25, 25, 45 





00059000 

25 DHE = DENHEI I ) 





00060000 

GO TO 50 





0006 1000 

45 DHE = DJHE ( ALT ) 





00062000 

50 RETURN 





00063000 

END 





00064000 

FUNCTION T ( ALT ) 





00065000 

IMPLICIT REAL*8( A-H.O-Z ) 





00066000 

REAL*8 





00067000 

1 MASS, MM, MN2 , M02 , HOI, MHE, MhAR 





00068000 

TEMPERATURE FROM bOKM TO TOP 





00069000 

COMMON/CALC/DR, CS, RE, ALTO, TINF, EDC, 

TLB, 

ZLB, 

MN2, 

MO 2, 

00070000 

1 MO 1 , MHE 





00071000 

COMMON /MOD/ T EM P ( 4 1 ) , MM(41), DENN2I41), 

DENU2 ( 41 ) , 

DEN01 (41 ) , 

00072000 

1 DENHE (41), DENA ( 4 1 ) 





00073000 

I = (ALT -ALTO) / DR + .5 





00074000 

I F ( ALT- 80.E05) 15,15,20 





00075000 

15 T = 186.0 





00076000 

GO TO 50 





00077000 
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20 IF (ALT - ZLB ) 25, 25, 45 
25 T = TEMPI I) 

GO T i) 50 
45 T = TJ(ALT) 

50 RETURN 
END 

FUNCTION D JN2 ( ALT ) 

IMPLICIT RE A L*8 ( A-H , 0— Z ) 

J65 N2 DENSITY 

COMMON /C ALC/DR , CS, KE, ALTO, TINF, EDC , TLB, ZLB, M N 2 , 
1 HOI , MHE 
REAL**8 

1 MASS, MM, M I'M 2 , M02 , I'lpl, MHE, MBAK 

GLB = 980.665 / ((1. + ZLB' / KE ) =P«=2 ) 

ZETA = (ALT - ZLB) * (RE + ZLB) / (RE + ALT) 

X = (TINF-800.) / (750. +-1.722E-04 * ( T INF-800. ) **2 ) 

A = 1.- TLB / TINF 
S = 0.0291 **DEXP ( - X * X / 2.) 

SIGMA = ( S + 1.50 E-4) * l.E-5 
EXPSZ =DEXP( -SIGMA * ZETA) 

GAMMA = MN2 * GLB / (SIGMA * 8.314E 07 * TINF) 

DJN2 = 4 . 008 E 11 * ((1. - A) / (1. - A * EXPSZ)) 4* (1. 
1 DEXP (-SIGMA * GAMMA * ZETA) 

RETURN 

END 

FUNCTION DJ02 ( ALT ) 

IMPLICIT REAL**H( A-H.O-Z ) 

J 6 5 02 DENSITY 

C OM MON /CALC/DR, CS, RE, ALTU, TINF, EDC, TLB, ZLB, MN2, 
1 M 01, MHE 
REALMS 

1 MASS, MM, MN2 , H02 , M01, MHE, MBAR 

GLB = 980.665 / ((1. + ZLB / RE)***2) 

ZETA = (ALT - ZLB) * (RE + ZLB) / (RE + ALT) 

X = (TINF-800.) / (750. + 1.722E-04 * (TINF-800. ) ***2 ) 

A = 1.- TLB / TINF 
S = 0.0291 * D E X P ( — X * X / 2.) 

SIGMA = ( S + 1.50 E-4) * l.E-5 
EXPSZ =DEXP( -SIGI'.A * ZETA) 

GAMMA = M02 4 GLB / (SIGMA * 8.314E 07 * TINF) 

DJ02 = 7.495E 10 * ( ( 1 . - A ) / (1. - A * EXPSZ)) ** (1. 
1 DEXP (-SIGMA * GAMMA * ZETA) 

RETURN 


00U78000 
00079000 
00080000 
00081000 
00082000 
00083000 
00084000 
00085000 
00086000 
H 0 2 , 00087000 

00088000 
00089000 
00090000 
00091000 
00092000 
00093000 
00094000 
00095000 
00096000 
00097000 
00098000 
+ GAMMA) **00099000 
00100000 
00101000 
00102000 
00103000 
00104000 
00105000 
MO 2 , 00106000 

00107000 
00108000 
00109000 
00110000 
00111000 
00112000 
00113000 
00114000 
00115000 
00116000 
00117000 
+ GAMMA) **00118000 
00119000 
00120000 


END 

FUNCTION DJOl(ALT) 

IMPLICIT REAL*8( A-H.O-Z ) 

J65 n DENSITY 

COMMON/CALC/DR, CS, KE, ALTO, TINF, EDC, TLB, ZLB, MN2 , M02, 
1 Mill, MHE 
REAL**8 

1 MASS, MM, MN2, M02, HOI, MHE , MBAR 

GLB = 980.665 / ((1. + ZLB / RE>**2) 

ZETA = (ALT - ZLB) =* (KE + ZLB) / (RE + ALT) 

X = (TINF-800.) / (750. + 1.722E-04 * (TINF-800. ) *P*s= 2 ) 

A = 1 . — TLB / TINF 
S = 0.0291 =*DEXP (-X =* X / 2.) 

SIGMA = ( S + 1.50 E-4).* l.E-5 
EXPSZ =DEXP(— SIGMA * ZETA) 

GAMMA = MO 1 * GLB / (SIGMA * 8.314E 07 ** TINF) 


00121000 

00122000 

00123000 

00124000 

00125000 

00126000 

00127000 

00128000 

00129000 

00130000 

00131000 

00132000 

00133000 

00134000 

00135000 

00136000 
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O.lni = 7.600F 10 * ((1. - A) / (1. - A ^ EXPSZ)) ** (1. + GAMMA ) ¥0013700(1 


1 DEXP (-SIGMA ¥ GAMMA * ZETA) 

RETURN 

END 

FUtMCT I ON DJHE(ALT) 

IMPLICIT R E A L*8 ( A— H *0- Z ) 

J65 HE DENSITY 

COMMON /CALC/ OR, CS, RE, ALTO, TINE, EDC , TLB, ZLB, MM2, HU 2 , 

1 ■ K 0 1 , M H E 
REAL¥8 

1 MASS, MM, MN2, M02, MOl, MHE , MBAR 

ALPHA = -0.4 

GLB = 980.665 / ((1. + ZLB / RE)**2) 

ZETA = (ALT - ZLB) ¥ (KE + ZLB) / (RE + ALT) 

X = (TINF-800.) / (750. + 1.722E-04 ¥ ( T I I'M F — WOO- ) **2 ) 

A = 1.- TLB / T I N F 
S = 0.0291 ¥ 0 E X P ( — X ¥ X / 2.) 

SIGMA = ( S + 1.50 E -A ) * l.E-5 
EXPSZ =DEXP( -SIGMA ¥ ZETA) 

GAMMA = MHE ¥ GLB / (SIGMA ¥ H.314E 07 ¥ TINE) 

DJHE = 2.400E 07 ¥ ( ( 1. - A ) / ( 1. - A ¥ EXPSZ)) ¥* (1. + ALPHA 
1 + GAMMA) ¥DEXP(-SIGMA ¥ GAMMA ¥ ZETA) 

RETURN 


(■0133000 
00139000 
00140000 
00141000 
00142000 
00143000 
00144000 
UO 14 5000 
00146000 
00 147000 
OU 1480 00 
00149000 
UO 150000 
00151000 
00152000 
00153000 
00154000 
00 155000 
00156000 
00157000 
00158000 
00159000 


END 

FUNCTION TJ (ALT) 

IMPLICIT REAL*8(A-H,0-Z) 

R E A L ¥ 8 

1 MASS* MM, MN2 , MO 2 , HOI, MHE, MBAk 

J65 TEMPERATURE 

COM MON /CALC/DR, CS, RE, ALTO, TINF, EDC, TLB, ZLB, MN2, M02, 
1 M01 , MHE 

X = (TINF-800.) / (750. + 1.722E-04 * ( T INF-800. )*#2 . ) 

S = 0.0291 ¥OEXP( -X¥X/2. ) 

SIGMA = ( S + 1.50 E-4 ) ¥ l.E-5 

ZETA = (ALT - ZLB) ¥ (RE + ZLB) / (RE + ALT) 

TJ = TINF - (TINF - TLB) ¥QEXP (-SIGMA ¥ ZETA) 

RETURN 

END 


00160000 

00161000 

00162000 

00163000 

00164000 

00165000 

00166000 

00167000 

00168000 

00169000 

00170000 

00171000 

00172000 

00173000 

00174000 
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